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Prefacio 


La necesidad de acudir a herramientas estadfsticas para el analisis de datos en todas 
las areas del conocimiento, ha hecho que aparezcan con el correr de los anos nuevas 
metodologfas que, no obstante se centran en fundamentos probabilfsticos comunes, son 
especfficas para cada una de las diversas disciplinas del saber. Algunos ejemplos son, entre 
otros, la econometrfa, psicometrfa o la bioestadfstica. La gran relevancia que tiene 
actualmente a nivel mundial el tema ambiental ha hecho que los profesionales en estadfstica 
encaminen esfuerzos en el desarrollo de nuevas tecnicas apropiadas para el analisis de 
informacion enmarcada dentro de este contexto. Como consecuencia de este impulso 
surgio una nueva rama de la estadfstica, denominada environmetrics (estadfstica 
ambiental). Dentro de esta ultima, los metodos geoestadfsticos juegan un papel 
preponderante. 

El presente documento tiene como proposito servir de consulta a geologos, 
biologos, ecologos, agronomos, ingenieros, meteorologos y todos aquellos profesionales 
que se encargan del estudio de informacion ambiental georreferenciada. Se toma como 
base para las aplicaciones informacion de variables fisicoqufmicas y biologicas medidas en 
un estuario ubicado en la costa norte de Colombia. La razon fundamental para lo anterior, 
es que este escrito es uno de los resultados centrales de un proyecto de investigacion 1 , cuyo 
objetivo fundamental fue el de evaluar la aplicabilidad de algunos procedimientos 
estadfsticos en el analisis de datos medidos en este tipo de ecosistemas. 

El documento tiene un enfoque teorico-practico. Para el seguimiento completo de la 
teorfa descrita se requiere tener conocimientos basicos de algebra de matrices y de 
estadfstica matematica. Sin embargo aquellas personas que esten poco familiarizadas con 
estos temas, podran obviar la lectura de algunas secciones en las que se hacen desarrollos 
teoricos y centrar su atencion en la filosoffa de los metodos presentados y en las 
aplicaciones mostradas en cada uno de los capftulos del documento. Una resumen no 
exhaustivo de conceptos de algebra lineal y de estadfstica es hecho al final en el apendice. 

No obstante en el escrito se cubren diversos temas geoestadfsticos y se hacen 
aplicaciones de metodos recientes, es necesario acudir a la lectura de artfculos cientfficos y 
textos avanzados para lograr un buen dominio de esta metodologfa. Un libro formal desde 
el punto de vista matematico con aplicaciones en diversas disciplinas es Cressie (1993). 
Otras referencias pueden ser tomadas de la bibliograffa. 


1 Proyecto "Analisis y aplicacion de tecnicas geoestadfsticas en la modelacion de procesos estocasticos 
relacionados con variables ecologicas en ambientes estuarinos", cofinanciado por INVEMAR y 


COLCIENCIAS. 
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Introduction 


El estudio de fenomenos con correlacion espacial, por medio de metodos 
geoestadfsticos, surgio a partir de los anos sesenta, especialmente con el proposito de 
predecir valores de las variables en sitios no muestreados. Como antecedentes suelen 
citarse trabajos de Sichel (1947; 1949) y Krige (1951). El primero observo la naturaleza 
asimetrica de la distribucion del contenido de oro en las minas surafricanas, la equiparo a 
una distribucion de probabilidad lognormal y desarrollo las formulas basicas para esta 
distribucion. Elio permitio una primera estimacion de las reservas, pero bajo el supuesto de 
que las mediciones eran independientes, en clara contradiccion con la experiencia de que 
existen “zonas” mas ricas que otras. Una primera aproximacion a la solucion de este 
problema fue dada por geologo G. Krige que propuso una variante del metodo de medias 
moviles, el cual puede considerarse como el equivalente al krigeado simple que, como se 
vera mas adelante, es uno de los metodos de estimacion lineal en el espacio con mayores 
cualidades teoricas. La formulacion rigurosa y la solucion al problema de prediccion 
(estimacion en muchos textos geoestadfsticos) vino de la mano de Matheron (1962) en la 
escuela de minas de Paris. En los anos sucesivos la teorfa se fue depurando, ampliando su 
campo de validez y reduciendo las hipotesis necesarias (Samper y Carrera, 1990). De la 
minerfa las tecnicas geoestadfsticas, se han "exportado" a muchos otros campos como 
hidrologfa, ffsica del suelo, ciencias de la tierra y mas recientemente al monitoreo 
ambiental y al procesamiento de imagenes de satelite. 

Aunque la aplicacion de la herramienta geoestadfstica es bastante reciente, son 
innumerables los ejemplos en los que se ha utilizado esta tecnica en estudios ambientales 
con el animo de predecir fenomenos espaciales (Robertson, 1987; Cressie y Majure, 1995; 
Diggle et al., 1995). La columna vertebral del analisis geoestadfstico es la determinacion 
de la estructura de autocorrelacion entre los datos y su uso en la prediccion a traves de las 
tecnicas conocidas como kriging y cokriging. Otros temas importantes dentro del estudio 
de informacion georreferenciada son el diseno de redes de muestreo (McBratney et al., 
1981), la geoestadfstica multivariada (Wackernagel, 1995) y la simulacion (Deutsh y 
Joumel, 1992). 

La geoestadfstica es solo una las areas del analisis de datos espaciales. Es 
importante reconocer cuando la informacion georreferenciada es susceptible de ser 
analizada por medio de dicha metodologfa. Por ello en el documento se hace inicialmente 
una definicion global de estadfstica espacial y se describen las caracterfsticas especiales 
que enmarcan cada una de sus areas. 

En el estudio de informacion georreferenciada, de forma analoga a como se procede 
en la aplicacion de muchos procedimientos estadfsticos, la primera etapa que se debe 
cumplir es la del analisis exploratorio de datos (AED). Esta busca identificar localizacion, 
variabilidad, forma y observaciones extremas. Por ello en el primer capftulo del escrito se 
hace una revision de metodos empleados en el AED y se describen algunos particularmente 
utiles en el contexto del analisis de informacion georreferenciada. Posteriormente en el 
segundo capftulo, entrando en materia, se hace definicion de conceptos basicos dentro de la 
teorfa geoestadfstica. 
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En el tercer capftulo se describen los procedimientos empleados para identificar de 
manera experimental (con base en datos muestrales) la estructura de autocorrelacion 
espacial, para algunas distancias dadas, de un conjunto de datos de una variable. Se muestra 
tambien como generalizar dicha estructura para cualquier distancia entre los sitios de 
observacion. Una vez detectada la autocorrelacion espacial, el siguiente paso es la 
prediccion en sitios de la region de estudio donde no se ha hecho medicion de la variable de 
interes. Esto es llevado a cabo por medio de alguno de los procedimientos kriging que son 
descritos en el capftulo cuatro. Por ultimo, en el capftulo cinco, se hace referenda a temas 
especiales dentro del analisis geoestadfstico como cokriging, componentes principales 
regionalizados, diseno de redes de muestreo y simulacion. En cada seccion del documento, 
despues de que han sido expuestos los aspectos teoricos esenciales de cada tecnica, se 
muestran aplicaciones practicas. 
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Capitulo Uno 

Datos Espaciales y Analisis Exploratorio 


En las secciones 1.1 y 1.2 se define estadfstica espacial y se mencionan sus 
subdivisiones. Lo anterior se hace con el proposito unico de que el lector identifique el 
alcance del tema considerado dentro del escrito. Por ello a partir de la seccion 1.3 de este 
capitulo y en los capftulos siguientes se consideran solo temas referentes a geoestadfstica 

1.1. Estadfstica Espacial. 

Estadfstica espacial es la reunion de un conjunto de metodologfas apropiadas para el 
analisis de datos que corresponden a la medicion de variables aleatorias en diversos sitios 
(puntos del espacio o agregaciones espaciales) de una region. De manera mas formal se 
puede decir que la estadfstica espacial trata con el analisis de realizaciones de un proceso 
estocastico{z(i) :«£/)}, en el que se R d representa una ubicacion en el espacio euclidiano d- 
dimensional, Z(s) es una variable aleatoria en la ubicacion s y s varfa sobre un conjunto de 
indices Dcz R d . 

1.2. Areas de la Estadfstica Espacial. 

La estadfstica espacial se subdivide en tres grandes areas. La pertinencia de cada una 
de ellas esta asociada a las caracterfsticas del conjunto D de indices del proceso estocastico 
de interes. A continuacion se mencionan dichas areas y se describen las propiedades de D 
en cada una de estas. 

Geoestadfstica: Las ubicaciones s provienen de un conjunto D continuo y son 
seleccionadas a juicio del investigador (D fijo). Algunos ejemplos de datos que pueden ser 
tratados con esta metodologfa son: Niveles de un contaminante en diferentes sitios de una 
parcela, contenidos aurfferos de una mina, valores de precipitacion en Colombia medida en 
las diferentes estaciones meteorologicas en un mes dado o los niveles piezometricos de un 
acuffero. En los ejemplos anteriores es claro que hay continuidad espacial, puesto que en 
cualquier sitio de la parcela, de la mina, de Colombia o del acuffero pueden ser medias las 
correspondientes variables. Es importante resaltar que en geoestadfstica el proposito 
esencial es la interpolacion y si no hay continuidad espacial pueden hacerse predicciones 
carentes de sentido. Por ejemplo si la variable medida es produccion de cafe en las fincas 
cafeteras del departamento del Quindfo, hacer interpolacion espacial y realizar un mapa de 
distribucion de la produccion cafetera puede ser carente de sentido porque podrfan hacerse 
predicciones sobre areas urbanas o no cultivadas con cafe. Ademas de lo anterior las 
mediciones, no obstante sean georreferenciadas, corresponden a una agregacion espacial 
(finca) mas que a un punto del espacio. En la parte de arriba, al comienzo de este parrafo, 
se menciono que D debfa ser fijo. A este respecto cabe aclarar que el investigador puede 
hacer seleccion de puntos del espacio a conveniencia o puede seleccionar los sitios bajo 
algun esquema de muestreo probabilfstico. 



• Lattices (enmallados): Las ubicaciones s pertenecen a un conjunto D discreto y son 
seleccionadas por el investigador (D fijo). Estas pueden estar regular o irregularmente 
espaciadas. Algunos ejemplos de datos en lattices son los siguientes: Tasa de morbilidad 
de hepatitis en Colombia medida por departamentos, tasa de accidentalidad en sitios de una 
ciudad, produccion de cana de azucar en el departamento del Valle del Cauca segun 
municipio, colores de los pixeles en interpretacion de imagenes de satelite. En los ejemplos 
anteriores se observa que el conjunto de ubicaciones de interes es discreto y que estas 
corresponden a agregaciones espaciales mas que a un conjunto de puntos del espacio. Es 
obvio que la interpolacion espacial puede ser carente de sentido con este tipo de datos. 

• Patrones Espaciales: las ubicaciones pertenecen a un conjunto D que puede ser 
discreto o continuo y su seleccion no depende del investigador (D aleatorio). Ejemplos de 
datos dentro de esta area son: Localizacion de nidos de pajaros en una region dada, puntos 
de imperfectos dentro de una placa metalica, ubicacion de los sitios de terremoto en 
Colombia o cuadrantes de una region con presencia de una especie particular. Debe notarse 
que en los ejemplos anteriores hay aleatoriedad en la seleccion de los sitios, puesto que la 
ubicacion de los nidos de los pajaros, de los imperfectos dentro de la placa metalica, de los 
sitios de terremoto o de los cuadrantes con presencia de la especie, no dependen del 
criterio del investigador. Una vez se ha hecho la seleccion de sitios es posible hacer 
medidas de variables aleatorias en cada uno de ellos. Por ejemplo si en primera instancia se 
establece la ubicacion de arboles de pino dentro de un bosque, es posible que sea de interes 
medir en cada uno de los arboles el diametro o la altura. En general el proposito de analisis 
en estos casos es el de determinar si la distribucion de los individuos dentro de la region es 
aleatoria, agregada o uniforme. 

1.3. Datos Georrferenciados 

Las mediciones de las caracterfsticas de interes en un estudio regionalizado tienen 
implfcitamente asociadas las coordenadas de los sitios en donde estas fueron tomadas. 
Cuando el area de estudio es considerablemente grande se usa un geoposicionador para 
establecer dichas coordenadas. En otros casos, por ejemplo en disenos experimentales con 
parcelas, es suficiente con hacer asignaciones segun pianos cartesianos. Un esquema 
general de datos georreferenciados es el siguiente: 


Sitio 

Latitud 

Norte 

Longitud 

Este 

X; 

x 2 



1 

- 

- 

Xu 

X]2 


%lp 

2 

- 

- 

X21 

X22 


%2p 

3 

- 

- 

X31 

X32 


%3p 

4 

- 

- 

X 4 1 

X42 


X4p 

n 

- 

- 

%nl 

%n2 


%np 


En la tabla anterior n es el numero de sitios muestreados y p el de variables medidas en 
cada uno de ellos. Cada xy corresponde a la medida de la variable X, (j = 1, 2,..., p ) en el 
sitio i ( i= 1, 2,..., n), que puede ser cuantitativa o categorica. Algunas de las variables 
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pueden estar mas intensamente muestreadas que las otras (Xy faltantes). Las coordenadas 
pueden ser planas, geograficas (grados, minutos y segundos) o cartesianas. Sin embargo la 
posible utilizacion de unas u otras depende del software empleado para los analisis. 

1.4. Justification del Analisis Exploratorio de Datos Espaciales . 

En la aplicacion de la geoestadfstica es de suma importancia, al igual que en otros 
procedimientos estadfsticos (por ejemplo los modelos ARIMA dentro de la teorfa de series 
de tiempo), el analisis grafico. La identificacion de valores extremos y su ubicacion 
geografica, la evaluacion de la forma de la distribucion y el calculo de medidas de 
localizacion, variabilidad y correlacion es muy importante para establecer si algunos 
supuestos necesarios para la aplicacion de la teorfa geoestadfstica son validos o para definir 
que procedimiento de prediccion es el mas conveniente. Por ejemplo, como se vera en el 
capftulo cuatro, la decision de usar kriging ordinario o kriging universal se fundamenta en 
identificar si la media es o no constante en la region. El uso de kriging log-normal se basa 
en un criterio empfrico relacionado con la forma asimetrica de la distribucion de los datos 
muestrales. La decision de emplear cokriging depende de la deteccion de asociaciones 
entre las variables. 

1.5. Graficos Exploratorios 

Al igual que en un estudio exploratorio clasico, cuando se dispone de informacion 
georreferenciada se pueden emplear histogramas, diagramas de tallos y hojas y de caja y 
bigotes (Hoaglin et al., 1983) con el proposito de identificar localizacion, variabilidad, 
forma y observaciones extremas. Adicionalmente los graficos de dispersion son muy utiles 
tanto para la deteccion de relaciones entre las variables como para la identificacion de 
tendencias en el valor promedio de la variable en la region (relacion entre la variable 
medida y las coordenadas geograficas). Un supuesto fundamental en el analisis 
geoestadfstico es que el fenomeno es estacionario, para lo cual, entre otros aspectos, el 
nivel promedio de la variable debe ser constante en todos los puntos del area de estudio. 
Una deteccion de tendencia en el grafico de dispersion puede ser una muestra de que no se 
satisface dicho supuesto. El grafico se construye tomando como eje de las abcisas la 
variable que representa la coordenada geografica y en el eje de las ordenadas la variable 
cuantitativa de estudio. La observacion de la nube de puntos resultante, incluso el ajuste de 
una lfnea de regresion, permite establecer de manera empfrica si existe dicha tendencia. Un 
grafico de dispersion entre valores de la variable separados por una distancia espacial dada 
(dispersograma rezagado ) es util en la deteccion de autocorrelacion espacial. Otro grafico 
que tradicionalmente se emplea en la descripcion de datos espaciales es el de datos 
clasificados segun puntos de referenda (media, mediana, cuartfles). Este permite comparar 
zonas del sistema de estudio respecto a las magnitudes de las variables. 
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1.6 Aplicacion: Estudio exploratorio de la distribution de datos fisicoquimicos y 
biologicos medidos en el estuario Cienaga Grande de Santa Marta en Marzo de 1997. 

Con informacion de las variables salinidad, seston (mg/1), nitritos (pmol/1), 
silicatos(pmol/l) y clorofila a (pg/1) medidas en una jomada de muestreo realizada en 
marzo de 1997 en el estuario Cienaga Grande de Santa Marta (CGSM)(Fig. 1), se realizo 
un estudio exploratorio de datos. Los resultados encontrados son descritos a continuacion: 
En primera instancia, se evidencia en el diagrama de caja (Fig. 2) y en el grafico de tallos 
y hojas (Fig. 3) que, con excepcion de la variable nitritos, existe un comportamiento 
simetrico en las distribuciones de los datos. Se observa tambien en estas figuras, que en 
todas las variables se presentan algunos valores “atlpicos” o muy alejados del 
comportamiento general antes mencionado. Fo anterior, antes de ser tornado como un 
indicador de alta variabilidad o de errores de medicion, puede ser considerado como un 
reflejo del comportamiento espacial de las variables dentro del ecosistema. Fa simetrfa de 
la mayorfa de las variables hace pensar que existe una gran zona en donde las condiciones 
del sistema respecto a la calidad del agua son bastante similares (esto podrfa ser lo que se 
conoce como cuerpo de agua de la CGSM) y los valores “alejados” pueden estar 
representando las condiciones de sitios especificos, particularmente especiales dentro del 
sistema, como son la zona mas estuarina (sitios de muestreo cercanos al sitio Boca de la 
Barra, Fig. 1) y las de desembocaduras de los rfos que bajan de la Sierra Nevada de Santa 
Marta (costado oriental y sur del sistema, Fig. 1). 


74 30 ' 



Figura 1 . Area de estudio y cuadrfculas en que fue subdividido el sistema Cienaga Grande de Santa Marta 
para realizar la toma de muestras. Cada una de las 115 cuadrfculas tiene un area de 4 km2. Los datos fueron 
tornados en el centra de cada una de ellas. 

Fa afirmacion de que no existen problemas de alta variabilidad y que por el 
contrario los datos medidos son bastante homogeneos, puede confirmarse con los valores 
de los coeficientes de variacion (tabla 1). En su mayorfa estos son menores del 30% y por 
consiguiente indicadores de poca heterogeneidad en la informacion. 
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Tabla 1. Medidas de localization y variabilidad de algunas variables medidas en la superficie de la columna 
de agua del estuario Cienaga Grande de Santa Marta en Marzo de 1997. 


Medida 

Salinidad 

Seston (mg/1) 

Nitritos (umol/1) 

Silicatos 

(umol/1) 

Clorofila a 
(ug/l) 

Media 

17.6 

218.28 

0.436 

244.94 

132.44 

Mediana 

16.9 

215 

0.350 

251.83 

137.37 

Minimo 

13.02 

103 

0.01 

10.98 

2.91 

Maximo 

34.9 

318 

1.61 

358 

198.3 

Cuartfl Inferior. 

15.97 

191 

0.210 

226.52 

124.43 

Cuartfl Superior. 

18.04 

248 

0.6 

278.43 

149.29 

Desviacion Estandar 

2.79 

41.1 

0.309 

61.43 

31.30 

Coeficiente de Variation 

16.1 

18.8 

70.8 

25.07 

23.7 


Las medidas de localizacion (media y mediana, tabla 1) toman valores similares a 
los reportados en otros estudios para la misma epoca del ano. Una discusion a este respecto 
se encuentra en Hernandez (1986) y Hernandez y Gocke (1990). 


7 
5 
3 
1 

-1 
-3 
-5 

Salinidad Seston Nitritos Silicatos Clorofila 

Variable 

Figura 2. Diagramas de caja de algunas variables medidas en la superficie de la columna de agua del estuario 
Cienaga Grande de Santa Marta en Marzo de 1997. Las variables fueron estandarizadas antes de construir los 
diagramas. 

El grafico de dispersion de la variable salinidad (una de las de mayor relevancia en 
el establecimiento del comportamiento espacial de las variables en el sistema) respecto a las 
coordenadas latitud y longitud (Fig. 4), permite apreciar una leve tendencia en la magnitud 
de la variable a lo largo de estas direcciones, lo que hace suponer que, a pesar de la 
homogeneidad antes mencionada, el valor promedio de la misma no es constante en toda la 
region. Lo anterior se puede comprobar en el grafico 5, en donde se aprecia que en una 
gran parte de la zona centra de la Cienaga y hacia la desembocadura de los rfos Sevilla y 
Aracataca la magnitud de la variable es menor a la de los restantes sitios de muestreo. Esta 
figura revela claramente la influencia que tienen las entradas de agua (tal vez exceptuando 
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la entrada del no Fundacion) en el comportamiento de esta variable. Los valores 
relativamente altos, respecto a los antes descritos, en la zona occidental pueden ser 
consecuencia del proceso de lavado de suelos hipersalinos que se da en epoca de lluvias en 
el complejo Pajarales (sistema con el que tiene frontera la Cienaga) y que llegan al sistema 
a traves de los Canales Grande y Cl arm (Fig. 5). Los valores "altos" en la zona sur pueden 
ser de igual forma causados por la influencia del canal Grande y por circulacion de las 
masas de agua dentro del sistema (contrario a las manecillas del reloj) 


13 



a). Salinidad 


b). Seston 


Bajo 

127 

2 13° 

5 

4 14* 

14 

11 14° 

6677889 

18 15* 

0011224 

30 15° 

566888899999 

46 16* 

0011222333333334 

(13) 16° 

5555666777899 

55 17* 

011112233344 

43 17° 

5556777788889 

30 18* 

0000122233344 

17 18° 

555556777888 

5 19* 

01 

Alto 

283,332,342 


c). Nitritos 

5 0 16888 

26 1 000022555777779999999 

43 2 11113333555558888 

(19) 3 0000022244444666699 

50 4 1133355555777777 

34 5 02248 

29 6 000033355599 

17 7 1668 

13 8 027 

10 9 366 

7 10 2 

6 11 3 

Alto 120,126,131,141,161 
e). Clorofila a. 

Baio 2.11,17,43,54,58,75 

8 9 7 

10 10 04 

21 11 01233334579 

41 12 00033444666677788899 

(18) 13 001124444556788889 

47 14 001111123334444567899 

26 15 0011333444467788 

10 16 0266 

6 17 00138 

Alto 198 


Bajol 10 
2 1 * 1 

3 IT 2 

5 IF 55 

15 IS 6667777777 

36 lo 888888888899999999999 

20 ) 2 * 00000001111111111111 

47 2T 2222222222222333333 

27 2F 44555555 

19 2S 6666777777777 

6 2o 99 

4 3* 0011 


d). Silicatos 

Bajo 1,1,1,1,2 
6 IS 6 

11 lo 88999 

22 2 * 00000111111 

43 2T 222222223333333333333 

(24 2F 444444444555555555555555 

45 2S 666666666666677777 

27 2o 888888889999 

15 3* 000000111 

6 3T 2223 

2 3F 4 

Alto 35 


Figura 3. Diagramas de tallos y hojas de algunas variables medidas en la superficie de la columna de agua 
del estuario Cienaga Grande de Santa Marta en Marzo de 1997. 
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Figura 4. Graficos de dispersion de valores de salinidad respecto a las coordenadas geograficas de medicion. 
Datos tornados en la superficie de la columna de agua del estuario Cienaga Grande de Santa Marta en Marzo 
de 1997. 
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• 13.00 to 16.40 
O 16.40 to 17.30 
□ 17.30 to 18.40 
^ 18.40 to 35.00 


Figura 5. Clasificacion de observaciones de la variable salinidad en intervalos (segun cuartfles) y ubicacion 
de estas dentro del area de estudio. Datos medidos en el estuario Cienaga Grande de Santa Marta en marzo de 


1997. 


Tabla 2. Matriz de correlacion calculada con base en informacion de algunas variables fisicoqui'micas y 
biologicas medidas en el estuario Cienaga Grande de Santa Marta en marzo de 1997. Los coeficientes que 
aparecen en negrita son significativos. 



Salinidad 

Seston 

Nitritos 

Silicatos 

Clorofila a 

Salinidad 

1 

-0.09 

-0.10 

-0.60 

-0.47 

Seston 


1 

-0.33 

0.06 

0.45 

Nitritos 



1 

-0.08 

-0.23 

Silicatos 




1 

0.46 

Clorofila a 





1 


Por ultimo de la matriz de correlacion (tabla 2) es posible afirmar que la abundancia 
fitoplanctonica, evaluda a traves de la concentracion de clorofila a, presenta correlacion 
significativa con las variables fisicoqufmicas medidas. Este patron de correlacion entre 
variables bioticas y abioticas en otros trabajos de menor intensidad muestral no ha podido 
ser detectado. En general los estudios realizados en la Cienaga Grande de Santa Marta en 
los que se pretende determinar los patrones de asociacion entre las variables biologicas y 
fisicoqufmicas siempre conducen a que la salinidad es la variable de mayor influencia en el 
regimen de productividad del sistema. Sin embargo estos resultados en primera instancia 
pueden estar detectando otro tipo de asociaciones. En los capftulos subsiguientes, cuando se 
realicen los mapas de distribucion espacial, se podran tener mas herramientas para discutir 
respecto a este tema. 
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Capitulo Dos 

Definiciones Basicas de Geoestadistica 


2.1. Definition de Geoestadistica 

La geoestadistica es una rama de la estadistica que trata fenomenos espaciales 
(Journel & Huijbregts, 1978). Su interes primordial es la estimacion, prediccion y 
simulacion de dichos fenomenos (Myers, 1987). Esta herramienta ofrece una manera de 
describir la continuidad espacial, que es un rasgo distintivo esencial de muchos fenomenos 
naturales, y proporciona adaptaciones de las tecnicas clasicas de regresion para tomar 
ventajas de esta continuidad (Isaaks & Srivastava, 1989). Petitgas (1996), la define como 
una aplicacion de la teorfa de probabilidades a la estimacion estadistica de valuables 
espaciales. 

La modelacion espacial es la adicion mas reciente a la literatura estadistica. 
Geologia, ciencias del suelo, agronomia, ingenieria forestal, astronomia, o cualquier 
disciplina que trabaja con datos colectados en diferentes locaciones espaciales necesita 
desaiTollar modelos que indiquen cuando hay dependencia entre las medidas de los 
diferentes sitios. Usualmente dicha modelacion concieme con la prediccion espacial, pero 
hay otras areas importantes como la simulacion y el diseno muestral (Cressie, 1989). 

Cuando el objetivo es hacer prediccion, la geoestadistica opera basicamente en dos 
etapas. La primera es el analisis estructural, en la cual se describe la correlacion entre 
puntos en el espacio. En la segunda fase se hace prediccion en sitios de la region no 
muestreados por medio de la tecnica kriging (capitulo 4). Este es un proceso que calcula un 
promedio ponderado de las observaciones muestrales. Los pesos asignados a los valores 
muestrales son apropiadamente determinados por la estructura espacial de correlacion 
establecida en la primera etapa y por la configuracion de muestreo (Petitgas, 1996). Los 
fundamentos basicos de estas etapas son presentados a continuacion. 

2.2. Variable Regionalizada. 

Una variable medida en el espacio de forma que presente una estructura de 
correlacion, se dice que es una valuable regionalizada. De manera mas formal se puede 
definir como un proceso estocastico con dominio contenido en un espacio euclidiano d- 
dimensional R d , {Z(x) : x e D c R d }. Si d = 2, Z (x) puede asociarse a una valuable medida 
en un punto x del piano (Diaz-Frances, 1993). En terminos practicos Z(x) puede verse como 
una medicion de una valuable aleatoria (p.ej. concentracion de un contaminante) en un 
punto v de una region de estudio. 

Recuerdese que un proceso estocastico es una coleccion de valuables aleatorias 
indexadas; esto es, pai'a cada * en el conjunto de indices D, Z(x) es una valuable aleatoria. 
En el caso de que las mediciones sean hechas en una superficie, entonces Z(x) puede 
interpretarse como la variable aleatoria asociada a ese punto del piano (x representa las 
coordenadas, planas o geograficas, y Z la variable en cada una de ellas). Estas valuables 
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aleatorias pueden representar la magnitud de una vaiiable ambiental medida en un conjunto 
de coordenadas de la region de estudio. 

2.3. Momentos de una Variable Regionalizada 

Sea {Z(x) : x e D c R d }el proceso estocastico que define la vaiiable regionalizada. 
Para cualquier n puntos xj, x?, x n , el vector aleatorio Z(x) = [z(x,), Z(x 2 ), • • • ,Z(x n )] r 

esta definido por su funcion de distribucion conjunta 

F[z 1 ,z 2 ,—,z n ]=P[z{x l )<z 1 ,Z{x 2 )< z 2 ,—,Z(x„)< z n \ 

Conocidas las densidades marginales univariadas y bivariadas se pueden establecer los 
siguientes valores esperados (momentos univariados y bivariados): 

• E(z(x ,)) = m(x i ) 






V(z(x,))=E[z(x,)-m(x ! )] 2 =crf 

C(z{x i ), z(xj )) = E[z{x i ) - m(x i )][z(xj ) - m[xj )J : Funcion de autocovarianza 


y{z(x i ), z(x. )) = —e[z(xi ) - z(x .)] : Funcion de semivarianza 

j 2 1 


2.4. Estacionariedad 

La variable regionalizada es estacionaria si su funcion de distribucion conjunta es 
invariante respecto a cualquier translacion del vector h, o lo que es lo mismo, la funcion de 

distribucion del vector aleatorio Z(x) = [z(x l ), Z(x 2 ), • • • ,Z(x n )] r es identica a la del vector 

Z(x)= [z(xj +h), Z(x 2 +h), ■••,Z(x n +h)Y para cualquier h. La teoria geoestadfstica se 
basa en los momentos arriba descritos y la hipotesis de estacionariedad puede definirse en 
terminos de estos: 


2.4.1 Estacionariedad de Segundo Orden 

Sea {Z(x): xgDc R d } una variable regionalizada definida en un dominio D contenido 
en R d (generalmente una variable medida en la superficie de una region) se dice que Z(x) es 
estacionario de segundo orden si cumple: 

a. E [ Z(x)] = m, ke R, Vx e D czR d . 

El valor esperado de la variable aleatoria es finito y constante para todo punto en el 
dominio. 


b. COV [ Z(x) , Z(x+h)\ = C(h) < oc 

Para toda pareja {z(x), Z(x + h )}la covarianza existe y es funcion unica del vector de 
separacion h. 

En la figura 6 se muestra el grafico de una variable regionalizada estacionaria. 
Exceptuando fluctuaciones aleatorias, el valor promedio de la variable no muestra una 
tendencia definida en alguna direccion. 
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La existencia de la covarianza implica que la varianza existe, es finita y no depende de 
h, es decir V(z(x l )) = C(0) = cr . Asf mismo la estacionariedad de segundo orden implica 
la siguiente relacion entre la funcion de semivarianza y la de autocovarianza: 

y(z(x + h), Z(x)) = y(h) = ^E[z(x + h) - m - Z(x)+ mf 

= ^\e(z(x+ h)—m)~ +E(z(x) — m) 2 - 2 E(z(x + h) — m)(z(x)—m)\ 

= \ (j2 + \ cj2 - E ^ x + / 0- m)(z(x)- 7 n)} 

= a 2 -C{h). 



Figura 6. Representation de una superficie interpolada para una variable regionalizada estacionaria 

2.4.2. Estacionariedad Debit o Intrinseca 

Existen algunos fenomenos ffsicos reales en los que la varianza no es finita. En estos 
casos se trabaja solo con la hipotesis que pide que los incrementos [Z(x+h) - Z(x)] sean 
estacionarios, esto es (Clark, 1979): 

a. Z(x) tiene esperanza finita y constante para todo punto en el dominio. Lo que implica que 
la esperanza de los incrementos es cero. 

e[z(x + Ii) - Z(x)\ = 0 

b. Para cualquier vector h, la varianza del incremento esta definida y es una funcion unica 
de la distancia. 

V [ Z(x+h) - Z(x)] = E[Z(x + h ) - Z{x)f = 2 y (h) 
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Es claro que si una variable regionalizada es estacionaria fuerte entonces tambien sera 
estacionaria debil. El concepto de estacionariedad es muy util en la modelacion de series 
temporales (Box & Jenkins, 1976). En este contexto es facil la identificacion, puesto que 
solo hay una direccion de variacion (el tiempo). En el campo espacial existen multiples 
direcciones y por lo tanto se debe asumir que en todas el fenomeno es estacionario. Cuando 
la esperanza de la variable no es la misma en todas las direcciones o cuando la covarianza o 
correlacion dependan del sentido en que se determinan, no habra estacionariedad. Si la 
correlacion entre los datos no depende de la direccion en la que esta se calcule se dice que 
el fenomeno es isotropico , en caso contrario se hablara de anisotropia. En Isaaks y 
Srivastava (1989) se definen los posibles tipos de anisotropia y se proponen algunas 
soluciones. Cressie (1993) discute cual debe ser el tratamiento en caso de que la media no 
sea constante. 

En casos practicos resulta compleja la identificacion de la estacionariedad. Suelen 
emplearse graficos de dispersion de la variable respecto a las coordenadas, de medias 
moviles y de valores clasificados segun puntos de referencia, con el proposito de identificar 
posibles tendencias de la variable en la region de estudio. L a isotropia es estudiada a traves 
del calculo de funciones de autocovarianza o de semivarianza muestrales (capitulo3) en 
varias direcciones. Si estas tienen formas considerablemente distintas puede no ser valido el 
supuesto de isotropia. Finalmente una variable regionalizada sera no estacionaria si su 
esperanza matematica no es constante, esto es si E[z(x)]= m(x). En la figura 7 se 
representa una variable regionalizada en la que existe tendencia en el valor promedio de la 
variable, lo cual es claro indicador de no estacionariedad. 



Figura 7. Representation de una superficie interpolada para una variable regionalizada no estacionaria 
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Capitulo Tres 

Correlation Espacial Muestral y Ajuste de Modelos 


3.1. Funciones de Correlation Espacial 

La primera etapa en el desarrollo de un analisis geoestadfstico es la determinacion 
de la dependencia espacial entre los datos medidos de una variable. Esta fase es tambien 
conocida como analisis estructural. Para llevarla a cabo, con base en la informacion 
muestral, se usan tres funciones: El semivariograma, el covariograma y el correlograma. A 
continuacion se hace una revision de los conceptos asociados a cada una de ellas y se 
describen sus bondades y limitaciones. 

3.1.1. Variograma y Semivariograma. 

Cuando se definio la estacionariedad debil en el capitulo anterior se menciono que se 
asurma que la varianza de los incrementos de la variable regionalizada era finita. A esta 
funcion denotada por 2y(h) se le denomina variograma. Utilizando la definition teorica de 
la varianza en terminos del valor esperado de una variable aleatoria, tenemos: 

2y(h) = v{z(x+h)-Z(x)) 

=e({Z(x + h)-Z(x)) 2 ] - ( E{Z(x + h)-Z(x))f 

K -v- J 

0 

=l((Z(x + /z)-Z(x)) 2 ) 


La mitad del variograma y(h), se conoce como la funcion de semivarianza y 
caracteriza las propiedades de dependencia espacial del proceso. Dada una realization del 
fenomeno, la funcion de semivarianza es estimada, por el metodo de momentos, a traves 
del semivariograma experimental, que se calcula mediante (Wackemagel, 1995): 

Y(Z(x + /z)-Z(x)} 2 

r(h) = — --- 

2 n 

donde Z(x) es el valor de la variable en un sitio x, Z(x+h) es otro valor muestral separado 
del anterior por una distancia h y n es el numero de parejas que se encuentran separadas por 
dicha distancia. La funcion de semivarianza se calcula para varias distancia h. En la 
practica, debido a irregularidad en el muestreo y por ende en las distancias entre los sitios, 
se toman intervalos de distancia {[o, h\ (h, 2h\ (2h, 3h\ •■■}y el semivariograma experimental 
corresponde a una distancia promedio entre parejas de sitios dentro de cada intervalo y no a 
una distancia h especffica. Obviamente el numero de parejas de puntos n dentro de los 
intervalos no es constante. 

Para interpretar el semivariograma experimental se parte del criterio de que a menor 
distancia entre los sitios mayor similitud o correlacion espacial entre las observaciones. Por 
ello en presencia de autocorrelacion se espera que para valores de h pequenos el 
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semivariograma experimental tenga magnitudes menores a las que este toma cuando las 
distancias h se incrementan. 

3.1.2. Covariograma y Correlograma. 

La funcion de covarianza muestral entre parejas de observaciones que se encuentran a 
una distancia h se calcula, empleando la formula clasica de la covarianza muestral, por: 

n 

^(Z(x + h)—m)(Z(x) — m) 

C(h) = COV(Z(x + h), Z(x)) = - 

n 

Y j [Z(x + h)-Z{x)) 

= - m 2 = C(h) 

n 


donde m representa el valor promedio en todo punto de la region de estudio y n es el 
mimero de parejas de puntos que se encuentran a una distancia h. En este caso es tambien 
valida la aclaracion respecto a las distancias dadas en el ultimo parrafo de la pagina 
anterior. 

Asumiendo que el fenomeno es estacionario y estimando la varianza de la variable 
regionalizada a traves de la varianza muestral, se tiene que el correlograma muestral esta 
dado por: 

... COV{Z(x + h),Z(x)) C(h ) C(h) 

r{h)= -- — 7 — =- 

s x+h -s x S* C( 0 ) 

Bajo el supuesto de estacionariedad cualquiera de las tres funciones de dependencia 
espacial mencionadas, es decir semivariograma, covariograma o correlograma, puede ser 
usada en la determinacion de la relacion espacial entre los datos. Sin embargo como se 
puede observar en las formulas, la unica que no requiere hacer estimacion de parametros es 
la funcion de semivarianza. Por esta razon, fundamentalmente, en la practica se emplea el 
semivariograma y no las otras dos funciones. 

A continuacion se presenta un ejemplo ilustrativo del calculo de la funcion de 
semivarianza experimental: Suponga que se tienen medidas sobre una variable hipotetica 
cuyos valores estan comprendidos entre 28 y 44 unidades y su configuracion en una region 
de estudio es como se presenta en el esquema de la siguiente pagina. Como se indica en la 
representacion, la distancia entre cada par de puntos contiguos es de 100 unidades. Luego si 
existe un punto faltante la distancia entre los dos valores ubicados a cada lado de este sera 
de 200 unidades. Veamos como calcular bajo esta situacion el semivariograma 
experimental. Por simplicidad se calcularan solo los semivariogramas en sentido 
(izquierda-derecha) e (inferior-superior), debido a que para obtener un semivariograma 
experimental en el que solo se tenga en cuenta la distancia y no la orientacion 
(.semivariograma omnidirectional), se requerirfa calcular la distancia euclidiana para un 
mimero considerablemente alto de parejas de puntos. 
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En primer lugar en sentido izquierda-derecha se encuentran todas las parejas de puntos 
que estan a una distancia de 100 unidades y se calcula el semivariograma como: 

f (100) = (38 - 37) 2 + (37 - 35) 2 + (29 - 30) 2 + ... + (37 - 36) 2 12* 36 = 1.458 
analogamente para la distancia de 200 unidades 

y (200) = (40 - 44) 2 + (40 - 40) 2 + (42 - 39) 2 + ... + (29 - 32) 2 12* 36 = 3.303 
Similarmente se procede para otras distancias y para el sentido inferior-superior. Los 
valores calculados de el semivariograma se muestran en la siguiente tabla. 

Tabla 3. Valores de la funcion de semivarianza experimental en dos direcciones para el conjunto de datos 
hipoteticos. 


Distancia 

Semivarianza 

Sentido Izquierda-Derecha 

Semivarianza 

Sentido Inferior-Superior 

100 

1.45 

5.34 

200 

3.30 

9.87 

300 

4.31 

18.88 

400 

6.69 

27.53 


A1 graficar los valores de la funcion de semivarianza experimental dados en la tabla 
anterior (Fig. 8 ) se observa que en sentido inferior-superior el semivariograma es mayor 
que en sentido izquierda-derecha, luego la conclusion mas relevante para este conjunto de 
datos es que la estructura de correlacion espacial no solo depende de la distancia entre los 
sitios, sino de su orientacion. En otras palabras el fenomeno podrfa ser anisotropico. 
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Figura 8. Funcion de semivarianza experimental en dos direcciones para un conjunto de datos hipoteticos. 

3.2. Modelos Teoricos de Semivarianza. 

Como se vera a partir del capftulo cuatro la solucion del problema de prediccion 
espacial kriging requiere del conocimiento de la estructura de autocorrelacion para 
cualquier posible distancia entre sitios dentro del area de estudio. En la presentacion del 
semivariograma experimental dada anteriormente se indico que este es calculado solo para 
algunas distancias promedios particulares. Por ello se hace necesario el ajuste de modelos 
que generalicen lo observado en el semivariograma experimental a cualquier distancia. 
Existen diversos modelos teoricos de semivarianza que pueden ajustarse al semivariograma 
experimental. En Samper y Carrera (1990) se presenta una discusion respecto a las 
caracterfsticas y condiciones que estos deben cumplir. En general dichos modelos pueden 
dividirse en no acotados (lineal, logaritmico, potencial) y acotados (esferico, exponencial, 
gaussiano) (Warrick et al., 1986). Los del segundo grupo garantizan que la covarianza de 
los incrementos es finita, por lo cual son ampliamente usados cuando hay evidencia de que 
presentan buen ajuste. Todos estos modelos tienen tres parametros comunes (Fig. 9) que 
son descritos a continuacion: 

■ Efecto Pepita 

Se denota por Co y representa una discontinuidad puntual del semivariograma en el 
origen (Fig. 9). Puede ser debido a errores de medicion en la variable o a la escala de la 
misma. En algunas ocasiones puede ser indicativo de que parte de la estructura espacial se 
concentra a distancias inferiores a las observadas. 

■ Meseta 

Es la cota superior del semivariograma. Tambien puede definirse como el limite del 
semivariograma cuando la distancia h tiende a infinito. La meseta puede ser o no finita. Los 
semivariogramas que tienen meseta finita cumplen con la hipotesis de estacionariedad 
fuerte; mientras que cuando ocurre lo contrario, el semivariograma define un fenomeno 
natural que cumple solo con la hipotesis intrfnseca. La meseta se denota por Ci o por (Co + 
Ci) cuando la pepita es diferente de cero. Si se interpreta la pepita como un error en las 
mediciones, esto explica porque se sugiere que en un modelo que explique bien la realidad, 
la pepita no debe representar mas del 50% de la meseta. Si el ruido espacial en las 
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mediciones explica en mayor proporcion la variabilidad que la correlacion del fenomeno, 
las predicciones que se obtengan pueden ser muy imprecisas. En la figura 9 se representa 
este parametro para el caso de uno de los modelos acotados. 



• SEMEXP 
-MODELO 


Figura 9. Comportamiento tfpico de un semivariograma acotado con una representacion de los parametros 
basicos. SEMEXP corresponde al semivariograma experimental y MODELO al ajuste de un modelo teorico. 


■ Rango 

En terminos practicos corresponde a la distancia a partir de la cual dos observaciones 
son independientes. El rango se interpreta como la zona de influencia. Existen algunos 
modelos de semivariograma en los que no existe una distancia finita para la cual dos 
observaciones sean independientes; por ello se llama rango efectivo a la distancia para la 
cual el semivariograma alcanza el 95% de la meseta. Entre mas pequeno sea el rango, mas 
cerca se esta del modelo de independencia espacial. El rango no siempre aparece de manera 
explfcita en la formula del semivariograma. En el caso del modelo esferico (3.2.1), el rango 
coincide con el parametro a, que se utilizara en las ecuaciones mas adelante. Sin embargo, 
en el modelo exponencial (3.2.2), el rango efectivo es a/3 y en el modelo gaussiano (3.2.3) 
es a/ a/3. 


3.2.1. Modelo Esferico 

Tiene un crecimiento rapido cerca al origen (Fig. 10), pero los incrementos marginales 
van decreciendo para distancias grandes, hasta que para distancias superiores al rango los 
incrementos son nulos. Su expresion matematica es la siguiente: 


y(i>)= 


C 0 + C] 


Co + C] 


( 

3 

(h\ 

1 l \ 

(h^ 

3 ^ 

V 

V 


' 2' 

yr/J 

/ 


h < a 
h > a 


En donde C/ representa la meseta, a el rango y h la distancia. 
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3.2.2. Modelo Exponencial 


Este modelo se aplica cuando la dependencia espacial tiene un crecimiento exponencial 
respecto a la distancia entre las observaciones. El valor del rango es igual a la distancia para 
la cual el semivariograma toma un valor igual al 95% de la meseta (Fig. 10). Este modelo 
es ampliamente usado. Su expresion matematica es la siguiente: 


/(h)- C 0 + C l 


1-exp 


3 h 


a 


W 


J) 


3.2.3. Modelo Gaussiano 


Al igual que en el modelo exponencial, la dependencia espacial se desvanece solo en 
una distancia que tiende a infinito. El principal distintivo de este modelo es su forma 
parabolica cerca al origen (Fig. 10). Su expresion matematica es: 


+ Cl 


f-/, 2 Y) 

1-exp 

— 

2 


V 

l « )) 



-Esferico 

.Exponencial 

-Gaussiano 


Figura 10. Comparacion de los modelos exponencial, esferico y Gaussiano. La tinea punteada vertical 
representa el rango en el caso del modelo esferico y el rango efectivo en el de los modelos exponencial y 
gaussiano. Este tiene un valor de 210, respecto a una escala simulada entre 0 y 300. El valor de la meseta es 
30 y el de la pepita 0. El 95% de la meseta es igual a 28.5. 


3.2.4. Modelo Monomicos. 

Corresponden a los modelos que no alcanzan la meseta (Fig. 11). Su uso puede ser 
delicado debido a que en algunos casos indican la presencia de no estacionariedad en 
alguna direccion. 
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0 < 0 < 2 


Su formula matematica es la siguiente: 

y(h)= kh 9 

Obviamente cuando el parametro 0 es igual a uno el modelo es lineal y k representa la 
pendiente de la ecuacion de regresion con intercepto cero. Graficamente se pueden 
representar asr: 


CO 

N 

C 

.CO 

i_ 

CO 

> 

E 

CD 

c n 



1 <Theta <2 
-o— Theta = 1 
— OcTheta < 1 


Distancia (h) 

Figura 11. Comportamiento tipico de los modelos de semivarianza monomicos. 


3.2.5. Modelo de Independencia (Pepita Puro). 

Es indicativo de carencia de correlacion espacial entre las observaciones de una 
variable (Fig. 12). Es comun sumar este modelo a otro modelo teorico de semivarianza, 
para obtener lo que se conoce como semivariograma anidado. Lo anterior se sustenta en 
una propiedad de los semivariogramas que dice que cualquier combinacion lineal de 
semivariogramas con coeficientes positivos es un semivariograma. Su expresion 
matematica es: 

f 0 h = 0 

u a > donde Co>0 

I L. Q ll > 0 

Su representacion grafica es la siguiente: 

Kh).. 

Co - 


1 - +h 

Figura 12. Modelo de semivarianza teorico para variables sin correlacion espacial. 
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La estimation de los parametros de los modelos teoricos descritos puede ser llevada a 
cabo, entre otros metodos, por maxima verosimilitud (Cressie, 1993) o regresion no lineal 
(Gotway, 1991). Algunos paquetes de computo geoestadfsticos como GS+ (Gamma Design 
Software, 1999) traen incorporados procedimientos iterativos como el de Gauss-Newton 
para llevar a cabo la estimation. Otros como GeoEAS (Englund y Sparks, 1988) solo 
permiten el ajuste a sentimiento por el metodo de ensayo y error. 

Como se menciono en la section 4.2. cuando la autocorrelation no es igual en todas las 
direcciones entonces se dice que hay anisotropfa. Esa puede ser geometrica o zonal. La 
primera se presenta cuando los semivariogramas calculados en varias direcciones tienen 
igual meseta pero varfan en el rango. En el segundo caso todos los semivariogramas 
direccionales tiene igual rango pero diferente meseta. Algunas transformaciones apropiadas 
para solucionar la anistropfa y hacer valida la construction de un semivariograma 
omnidireccional se pueden encontrar en Isaaks y Srivastava (1989), Samper y Carrera 
(1990) y Cressie (1993). 

3.3. Aplicacion: Estimation de Modelos de Semivarianza para algunas variables 
fisicoquimicas y biologicas medidas en el estuario Cienaga Grande de Santa Marta. 

En esta section se hace una interpretation practica de resultados encontrados al hacer 
estimation de modelos teoricos de semivarianza para un conjunto de variables medidas en 
el estuario Cienaga Grande de Santa Marta (IGAC, 1973). Se consideran para el analisis 
datos tornados en dos niveles de la columna de agua (superficie y fondo), de las variables 
salinidad, oxfgeno disuelto (mg/1), solidos en suspension (mg/1), nitritos (pmol/lj y clorofila 
“a”(pg/lj, Ademas se estudian valores de profundidad (m) y transparencia (m). Para cada 
variable se obtuvieron 115 observaciones tomadas a lo largo de todo el sistema mediante 
muestreo sistematico de cuadrfculas (cuadrfculas de 4 km 2 )(Fig. 1). La ubicacion en cada 
punto de muestreo fue realizada mediante un geoposicionador GPS 100 SRVY II (Garmin, 
1993). Los semivariogramas experimentales y los correspondientes modelos estimados se 
hallaron por medio del software GS+ (Gamma Design Software. 1999). 

• Resultados y Discusion. 

Los semivariogramas experimentales encontrados (Figs. 13 y 14) indican que las 
variables presentan estructuras de dependencia espacial, puesto que en ningun caso la 
semivarianza es constante en funcion de la distancia. Los rangos encontrados en los 
modelos teoricos ajustados a los semivariogramas (tabla 4), superan los 11 km y en 
algunos casos este parametro alcanza los 25 km, lo cual resulta relativamente alto, teniendo 
en cuenta que la distancia entre los extremos sur y nolle del sistema (la mas amplia) no 
supera los 30 km. Lo anterior es un indicador de fuerte dependencia espacial para el caso 
considerado. Esto es sin duda conveniente puesto que desde un punto de vista teorico es 
conocido que un alto valor en el rango permite obtener curvas de prediction mas 
suavizadas reduciendo las magnitudes en varianzas de prediction (Dfaz-Frances, 1993) . Es 
importante resaltar respecto a los otros dos parametros, que en ningun caso el valor de la 
pepita supera el 50% del valor de la meseta (tabla 4), lo cual, segun Dfaz-Frances (1993), es 
recomendable para que el modelo de correlation espacial describa bien la realidad. Si el 
ruido espacial en las mediciones explica en mayor proportion la variabilidad que la 
correlation del fenomeno, las predicciones pueden ser muy imprecisas. 


28 




Distancia (m) 


Distancia (m) 



Distancia (m) 


Distancia (m) 



Figura 13. semivariogramas experimentales (calculados con los datos muestrales) y ajustes de modelos teoricos para las variables medidas en la superficie de la 
columna de agua de la Cienaga Grande de Santa Marta en marzo de 1997. a) salinidad; b)oxfgeno; c) solidos en suspension; d) nitritos; e) clorofila a; 
f)profundidad. 
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Figura 14. Semivariogramas experimentales (calculados con los datos muestrales) y ajustes de modelos teoricos para las variables medidas en el fondo de la 
columna de agua de la Cienaga Grande de Santa Marta en marzo de 1997. a) salinidad; b)oxfgeno; c) solidos en suspension; d) nitritos; e) clorofila a; f) 
transparencia. 
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Tabla 4. Modelos teoricos ajustados a semivariogramas experimentales de variables fisicoquimicas y 
biologicas medidas en dos niveles de la columna de agua de la Cienaga Grande de Santa Marta, durante una 
jornada de muestreo realizada en marzo de 1997. 


Variable 

Nivel 

Modelo 

Pepita 

Meseta 

Rango (in) 

7 

Salinidad 

Superficie 

Lineal 

0.179 

12.309 

20000 

0.89 


Fondo 

Lineal 

0.627 

11.752 

20000 

0.99 

Oxigeno 

Superficie 

Gaussiano 

1.830 

14.320 

12940 

0.99 


Fondo 

Esferico 

1.080 

4.211 

11650 

0.99 

Solidos en suspension 

Superficie 

Lineal 

1087 

1138 

22000 

0.90 


Fondo 

Lineal 

1408 

557 

20000 

0.67 

Nitritos 

Superficie 

Lineal 

0.071 

0.043 

22000 

0.87 


Fondo 

Lineal 

0.073 

0.077 

20000 

0.70 

Clorofila a 

Superficie 

Lineal 

389.2 

623.2 

18000 

0.91 


Fondo 

Lineal 

710 

616.4 

18000 

0.91 

Profundidad 


Gaussiano 

0.073 

0.121 

24850 

0.99 

Transparencia 


Gaussiano 

0.0069 

0.0019 

25000 

0.85 


Se puede afirmar que las variables oxigeno disuelto, profundidad y transparencia 
cumplen con la hipotesis de estacionariedad fuerte, dado que sus modelos son acotados 
(Biau et al., 1997; Samper y Carrera, 1990). De otro lado salinidad, solidos en suspension, 
nitritos y clorofila “a”, solo cumplen la hipotesis intrfnseca (estacionariedad debit) puesto 
que sus modelos son lineales (Evangelos y Flatman, 1988; Samper y Carrera, 1990). 

Debido a que los resultados arriba descritos respecto a los semivariogramas 
experimentales y al ajuste de modelos teoricos, confirman la hipotesis de autocorrelacion 
espacial en las caracterfsticas medidas en el estuario de estudio, es posible afirmar que los 
metodos geoestadfsticos pueden ser una herramienta de gran utilidad en la modelacion e 
interpretacion de fenomenos observados en este tipo de ecosistemas. Cuando se utilicen 
metodos estadfsticos tradicionales (regresion, analisis de varianza, tecnicas multivariadas, 
muestreo) para el analisis de este tipo de informacion, debe involucrarse en los 
correspondientes modelos la estructura de correlacion espacial implicita en los datos. 
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Capitulo Cuatro 

Prediction Espacial 


4.1. Prediction Espacial Optima. 

De la teorfa de la decision se conoce que si Z 0 es una cantidad aleatoria y Z'o es su 
predictor 2 , entonces L( Z 0 ;Z* 0 ) representa la perdida en que se incuiTe cuando se predice 

Z 0 con Z* 0 y el mejor predictor sera el que minimice e{l(z o ;Z^)/ z} con 
/ 

Z = {Zj, Z 2 , ■ ■ ■, Z n } , es decir el predictor optimo es el que minimice la esperanza 

conditional de la funcion de perdida. Si l{z {) ;Z 0 ) = [z 0 -Z* 0 }~ ^Z* 0 =e(z 0 iz). La expresion 
anterior indica que para encontrar el predictor optimo se requiere conocer la distribucion 
conjunta de la n+1 variables aleatorias. 

4.2. Definicion de Kriging. 

La palabra kriging 3 (expresion anglosajona) procede del nombre del geologo 
sudafricano D. G. Kiige, cuyos trabajos en la prediccion de reservas de oro, realizados en la 
decada del cincuenta, suelen considerarse como pioneros en los metodos de interpolacion 
espacial. Kriging encierra un conjunto de metodos de prediccion espacial que se 
fundamentan en la minimizacion del error cuadratico medio de prediccion. En la tabla 5 se 
mencionan los tipos de kriging y algunas de sus propiedades. En la secciones 4.3 y 4.4, se 
hace una presentacion detallada de ellos. 


Tabla 5. Tipos de predictores kriging y sus propiedades. 


TIPO DE 
PREDICTOR 

NOMBRE 

PROPIEDADES 

LINEAL 

• Simple 

• Ordinario 

• Universal 

• Son optimos si hay normalidad 
multivariada. 

• Independiente de la distribucion son los 
mejores predictores linealmente 
insesgados. 

NO LINEAL 

• Indicador 

• Probabilfstico 

• Log Normal, Trans- 

Gaussiano 

• Disyuntivo 

• Son predictores optimos. 


2 La palabra estimation es utilizada exclusivamente para inferir sobre parametros fijos pero desconocidos; 
prediction es reservada para inferencia sobre cantidades aleatorias. 

3 Algunos textos indican que en espanol la palabra adecuada serfa krigeado. 
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Los metodos kriging se aplican con frecuencia con el proposito de prediccion, sin 
embargo estas metodologfas tienen diversas aplicaciones, dentro de las cuales se destacan 
la simulation y el diseno de redes optimas de muestreo (capftulo 5). 

4.3. Kriging Ordinario 

Suponga que se hacen mediciones de la variable de interes Z en los puntos x h i= 1, 
2,..., n, de la region de estudio, es decir se tienen realizaciones de las variables Z(xj), . . . , 
Z(x n ), y se desea predecir Z(x 0 ), en el punto xo donde no hubo medicion. En esta 
circunstancia, el metodo kriging ordinario propone que el valor de la variable puede 
predecirse como una combination lineal de las n variables aleatorias asi: 

Z (xo) = Z] Z(xj) + A ,2 Z(x 2 ) + Z? Z(xj) + Z 4 Z(x 4 ) + Z 5 Z(xs) + . . . + Z„ Z(x n ) 

= XA-ZfXi) 

i=l 

en donde los ?q representan los pesos o ponderaciones de los valores originales. Dichos 
pesos se calculan en funcion de la distancia entre los puntos muestreados y el punto donde 
se va a hacer la correspondiente prediccion. La suma de los pesos debe ser igual a uno para 
que la esperanza del predictor sea igual a la esperanza de la valuable. Esto ultimo se conoce 
como el requisito de insesgamiento. 

Estadfsticamente la propiedad de insesgamiento se expresa a traves de: 

E(z‘(xJ=E(z(x 0 )) 


Asumiendo que el proceso es estacionario de media m (desconocida) y utilizando las 
propiedades del valor esperado, se demuestra que la suma de las ponderaciones debe ser 
igual a uno: 




V i=l 


= m 


^Z i E(z(x i )) = m 

i—l 



n n 

Zz=i 

!=1 !=1 

Se dice que Z*(xo) es el mejor predictor, lineal en este caso, porque los pesos se obtienen 
de tal manera que minimicen la vaiianza del error de prediccion, es decir que minimicen la 
expresion: 

v[z*{x 0 )-z{x 0 j) 


Aj = m => 
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Esta ultima es la caractenstica distintiva de los metodos kriging, ya que existen otros 
metodos de interpolacion como el de distancias inversas o el poligonal, que no garantizan 
varianza minima de prediccion (Samper y Carrera, 1990). La estimacion de los pesos se 

n 

obtiene minimizando V[Z*(x 0 )-Z(x 0 )] sujeto a ^ A, =1. 

i=l 


Se tiene que V[Z*(x 0 ) - Z(jc 0 )] = V[Z*(x 0 )] - 2COV[Z* (x 0 ),Z(x 0 )] + V[Z{x 0 )] 
Desagregando las componentes de la ecuacion anterior se obtiene los siguiente: 


V[Z (x 0 )]=V 


X^Zfx,) =J J Y J A i A j COV[Z(x i ),Z( Xj )] 

J=1 J i=l j=l 


En adelante se usara la siguiente notacion: COV\Z(x l ),Z(x j )] = C tj y V[ Z( x 0 )]=a 2 


De lo anterior COV [ Z*(x 0 ),Z(x 0 )]=COV 


2^AiZ(x,),Z( x 0 ) 

i=l 


n n 

=X ^COVlZix ,), Z(x 0 )] = X A,C l0 

i =1 i=1 


Entonces reemplazando, se tiene que: 


V[Z*(x 0 )-Z(x 0 )/=££A,AjG ij - 2 £M)v+<r 2 (0) 

i=l j=l i=l 

n 

Luego se debe minimizar la funcion anterior sujeta a la restriccion ^ A t =l . Este problema 

i=1 

de minimizacion con restricciones se resuelve mediante el metodo de multiplicadores de 
Lagrange. 


<4 =Z I Wr 2 X 4 c w+° 2+ 2 Ji 

i = l J-l * = 1 Multiplicador^_ 

de Lagrange 


f n 'N 

Zv/ 

i=l 


0 


Siguiendo el procedimiento acostumbrado para obtener valores extremos de una funcion, se 
deriva e iguala a cero, en este caso con respecto a A, y //: 


dA x 


f n \ 


(A} Cjj +2Aj Zv.,+ZZw#>- 2 Z ^c»+^ 2 + 2 H Z a '-' 


1=2 


i=2 j =1 


1=1 


V i=1 J 


dA x 


2A x C ll +2^AjC lj 

j=2 J 

-V-' 

0 


-2C l0 +2/i 


n n 

= 2y, A jC x j -2C 10 +2 n = 0 => Z ^jCij+M 

j=i i=i 


C l0 (1) 


34 



De manera analoga se determinan las derivadas con respecto a A 2 A,,: 

T/ 2\ n n 

= 2£AjC 2j -2C 20 +2jU = 0^ A J C 2J+M = C 20 (2) 

2 7=1 7=1 


d(oi ) 

dA„ 


2 Z A j C nJ- 2C »o +2 <“ = 0 => HtjCv+M = C «0 (3) 


7=1 7=1 

por ultimo derivamos con respecto a ju: 


dn 


n n 

2 ^ 2,-2 = 0 1(4) 


i=l 


i=l 


De (1), (2), (3), (4) resulta un sistema de (n + 1) ecuaciones con (n + 1) incognitas, que 
matricialmente puede ser escrito como: 


'C n ■ 

• c u r 

v 


(c 3 

WO 

C n i 

1 1 • 

C (m 1 

. . 1 oj 

K 


C„o 

l 1 J 


Ci, 


A = Cio 


por lo cual los pesos que minimizan el error de prediccion se determinan mediante la 
funcion de covariograma a traves de: 


A 


Cu 


Cio- 


Encontrando los pesos se calcula la prediccion en el punto x a . De foima analoga se procede 
para cada punto donde se quiera hacer prediccion. 

• Varianza de Prediccion del Kriging Ordinario 

Multiplicando (1), (2) y (3) por 2,- se obtiene: 


2 ; 


YjACij +M 

\J =1 


= 2 ; C a Vi, i = 1, 2, • • •, n. 


Sumando las n ecuaciones 


n n 


2 ,// = ^ 2 ,Q 

i =1 7=1 


i=l 


i=l 


i =1 7=1 i=l /=1 

Sustituyendo la expresion anterior en (0) 
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=° 2 

*=1 *=1 i= 1 

n 

=° 2 ~Yj X ‘ C ‘0~M ( 5 ) 

1=1 

• Estimation de Ponderaciones por medio de la Funcion de Semivarianza 

Los pesos A pueden ser estimados a traves de la funcion de semivarianza, para lo cual se 
requiere conocer la relacion entre las funciones de covariograma y de semivarianza. Antes 
de esto conveniente tener en cuenta la siguiente notacion: 

o 2 = V(Z(x)), Ytj = y(h), donde h es la distancia entre los puntos i y j y analogamente 
Cij = C(h). 

La relacion entre las dos funciones en cuestion es la siguiente: 

r^^Zix-) - Z(x,)) 2 } 

=^e[(Z(Xj)) 2 -2 (Z(xj)Z( Xi ) + (Z( X j )) 2 ] 

=^e[(Z( X j ) 2 ] - e[z(Xj)Z ( X[ )] + i e\z ( Xi )) 2 ] 

2 2 

= 4 l E )) 2 - ^ 2 1 + ^ x i )) 2 - ^ 2 1 - fc ( z ( L ' ) Z (- L )) - * ■ 2 1 

2 2 

=Uv(Z(x))\ +2 [v(Z(x))] - cov[z( Xj )Z ( X[ )] 

2 2 

= V[Z(X)\-C0V[Z(X ] )Z(X 1 )\ 

=a 2 -Cy ^ c y=a 2 - ru (6) 

Reemplazando (6) en (1), (2) y (3) se deteiTninan los pesos optimos X en terminos de la 
funcion de semivarianza: 

A/C,j + ^~ c io=i (v 2 -Yu)+v - (<7 2 -n o) 

1 7=1 7=1 

1 n 

=a2 Yj A j + 1* - cr2+ ^o 

7=1 7=1 

n n 

=cr 2 AjY\j +M-° 2 +Yi o => X A iYij ~ V =Yio 

7=1 7=1 
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Similarmente, 


dX 2 


X^r 2J -m = y 20 

j =i 


—^ = L X jrnj -M = 7n 0 


El sistema de ecuaciones se completa con (4). De acuerdo con lo anterior los pesos se 
obtienen en terminos del semivariograma a traves del sistema de ecuaciones: 



Para establecer la expresion de la correspondiente varianza del error de prediccion en 
terminos de la funcion de semivarianza se reemplaza (6) en (5), de donde: 

=cr 2 - YAj(a 2 - Yy ) +ju 

_ i =1 

a 2 k =a 2 -a^Ai+^ZiTij +M 

i =1 i=l 

n 

i =1 

Los pesos de kriging ordinario tambien pueden ser estimados mediante el uso del 
correlograma aplicando la siguiente relacion: p l] = C lJ /a 2 , caso en el que la correspondiente 
vaiianza de prediccion estaria dada por (Isaaks y Srivastava, 1989): 

f \ 

a t =(j2 

v J 

• Validacion del kriging. 

Existen diferentes metodos para evaluar la bondad de ajuste del modelo de 
semivariograma elegido con respecto a los datos muestrales y por ende de las predicciones 
hechas con kriging. El mas empleado es el de validacion cruzada, que consiste en excluir la 
observacion de uno de los n puntos muestrales y con los n-1 valores restantes y el modelo 
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de semivariograma escogido, predecir via kriging el valor de la variable en estudio en la 
ubicacion del punto que se excluyo. Se piensa que si el modelo de semivarianza elegido 
describe bien la estructura de autocorrelacion espacial, entonces la diferencia entre el valor 
observado y el valor predicho debe ser pequena. Este procedimiento se realiza en forma 
secuencial con cada uno de los puntos muestrales y asf se obtiene un conjunto de n “errores 
de prediccion”. Lo usual es calcular medidas que involucren a estos errores de prediccion 
para diferentes modelos de semivarianza y seleccionar aquel que optimice algun criterio 
como por ejemplo el del mfnirno error cuadratico medio (MECM). Este procedimiento es 
similar a la conocida tecnica de remuestreo Jacknife (Efron, 1982) empleada en diversos 
contextos estadfsticos para calcular varianzas de estimation, entre otros aspectos. Una 
forma descriptiva de hacer la validation crazada es mediante un grafico de dispersion de 
los valores observados contra los valores predichos. En la medida en que la nube de puntos 
se ajuste mas a una lfnea recta que pase por el origen, mejor sera el modelo de 
semivariograma utilizado para realizar el kriging. 

• Representation de las predicciones 

Una vez se ha hecho la prediction en un conjunto de puntos diferentes de los 
muestrales via kriging, se debe elaborar un mapa que de una representacion global del 
comportamiento de la variable de interes en la zona estudiada. Los mas empleados son los 
mapas de contornos, los mapas de residuos y los graficos tridimensionales. En el caso de 
los mapas de contornos, en primer lugar se divide el area de estudio en un enmallado y se 
hace la prediccion en cada uno de los nodos de este mismo. Posteriormente se unen los 
valores predichos con igual valor, generando asf las lfneas de contomo (isolfneas de 
distribucion). Este grafico permite identificar la magnitud de la variable en toda el area de 
estudio. Es conveniente acompanar el mapa de interpolaciones de la variable con los 
correspondientes mapas de isolfneas de los errores y de las varianzas de prediccion 
(posiblemente estimados a traves de metodos matematicos), con el proposito de identificar 
zonas de mayor incertidumbre respecto a las predicciones. 

• Intervalos de Confianza. 


Asumiendo que los errores de prediccion siguen una distribucion normal estandar y que 
son independientes, un intervalo de confianza del 100(l-a)%, 0<a <1, para Z(x) es: 


2 * M - Z^a/:>Z* (x) + Zj „/ a k 

percentil de una normal estandar. 


con z*(x) el valor calculado de la prediccion y z^ a/2 ti 
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• Ilustracion 

Suponga que se tiene una configuration de datos como la que se presenta en el 
esquema de abajo. Con base en siete datos observados (valores al lado del signo + por fuera 
de los tirculos numerados de 1 a 7) se quiere predecir un valor de la variable en el punto 
donde se encuentra el signo de interrogation, por fuera del circulo con el numero cero. 



La matriz de distancia euclidianas entre los sitios es la siguiente: 


litio 

0 

1 

2 

3 

4 

5 

6 

7 

0 

0.00 

4.47 

3.61 

8.06 

4.49 

6.71 

8.94 

13.45 

1 


0.00 

2.24 

10.44 

13.04 

10.05 

12.17 

17.80 

2 



0.00 

11.05 

13.00 

8.00 

10.05 

16.97 

3 




0.00 

4.12 

13.04 

15.00 

11.05 

4 





0.00 

12.37 

13.93 

7.00 

5 






0.00 

2.24 

12.65 

6 







0.00 

13.15 

7 








0.00 


Asumiendo que la estructura de correlation espacial entre los datos es estimada por un 
modelo exponencial y{h)= 10 (l - exp(-3/i/10))(pepita cero, meseta 10 y rango 10) o en 
terminos de la funcion de autocovarianza por c(h) = 10 (exp(- 3h/\ 0)), se encuentran las 
siguientes matrices que permiten encontrar los pesos para la prediction: 
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C„ C l2 C 13 Cj 4 c 15 c l6 c 17 

£■21 C 22 C 23 C 2 4 C 25 C 26 C 2 7 

C31 ^32 ^33 ^34 ^35 ^36 ^37 

C41 c 42 c 43 c 44 c 45 c 46 c 47 

Oil C 5 2 C 53 ^54 Q 5 ^56 ^57 

^61 c 62 C 63 C 64 C 65 C 66 C 67 

^71 C 72 ^73 C 74 C 7 5 C 76 C 77 

^1111111 


A f 10 5.11 0.44 0.20 0.49 0.26 0.05 

1 5.11 10 0.36 0.20 0.91 0.49 0.06 

1 0.44 0.36 10 2.90 0.20 0.11 0.36 

1 _ 0.20 0.20 2.90 10 0.24 0.15 1.22 

, “ 0.49 0.91 0.20 0.24 10 5.11 0.22 

1 0.26 0.49 0.11 0.15 5.11 10 0.19 

1 0.05 0.06 0.36 1.22 0.22 0.19 10 

oj [ 1 1 1 1 1 1 1 


1 

1 

1 

1 

1 

1 

1 

0 


0.127 -0.077 -0.013 -0.009 -0.008 -0.009 -0.012 0.136 

-0.077 0.129 -0.010 -0.008 -0.015 -0.008 -0.011 0.121 

-0.013 -0.010 0.098 -0.042 -0.010 -0.010 -0.014 0.156 

0.009 -0.008 -0.042 0.102 -0.009 -0.009 -0.024 0.139 

-0.008 -0.015 -0.010 -0.009 0.130 -0.077 -0.012 0.118 

-0.009 -0.008 -0.010 -0.009 -0.077 0.126 -0.013 0.141 

-0.012 -0.011 -0.014 -0.024 -0.012 -0.013 0.085 0.188 

0.136 0.121 0.156 0.139 0.118 0.141 0.188 -2.180 



con base en el vector estimado de parametros se encuentra que 

7 

Z 0 * =2>Z,. = (0.173)(477) + (0.318X696) + ... + (0.086)(0.18) = 592. 

1=1 

7 

con a\ =(7 2 -^Z,C !0 -jU = 10-[(0.173)(2.6l) + --- + (0.086)(0.18)]-0.907 
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4.4. Otros Metodos Kriging 

A continuation se mencionan algunos aspectos generates de otros metodos de 
prediccion espacial. Un estudio riguroso de ellos puede hacerse en Cressie (1993), Deutsch 
y Journel (1998) y Samper y Carrera (1990). 

4.4.1. Kriging Simple 

Suponga que hay una valuable regionalizada estacionaria con media (m) y covarianza 
conocidas. De manera analoga a como se define en modelos lineales (por ejemplo en 
diseno de experimentos) el modelo establecido en este caso es igual a la media mas un error 
aleatorio con media cero. La diferencia es que en este caso los errores no son 
independientes. 

Sea Z(x) la variable de interes medida en el sitio x. 

L'fz^)] = m 

Z (x) = m + e(x), con E^x)] = 0. 

El predictor de la valuable de interes en un sitio xo donde no se tiene information se define 
como: 

Z*(x 0 ) = m + £*{x 0 ), 

con £ (x 0 ) que corresponde a la prediccion del error aleatorio en el sitio xo. Despejando de 
la ecuacion anterior £ (x 0 ) = Z (x 0 ) - m. 

El predictor del error aleatorio se define por: 

£* Mi*, )=XA (Z (x, )-m). 

i =1 i=l 

de donde el predictor de la variable de estudio es: 

n n 

Z (x 0 )=m+ XZ,( z (a)- m ) = m + ) 

_ i =1 _ i =1 

El predictor es insesgado si: 

e[z * (x 0 )) = E(z(x 0 )) = m . Luego el predictor sera insesgado cuando e{£ (x 0 )) = 0. 

( \ n H 

E^(x 0 ))=X z ,4a) = X Z, (0) = 0. Por consiguiente, a diferencia del kriging ordinario, 

i =1 i =1 

en este caso no existen restricciones para las ponderaciones tendientes al cumplimiento de 
la condicion de insesgamiento. 
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La estimation de los pesos del metodo kriging ordinario se obtiene de tal forma que 
se minimice vV(*o)- £ (*o))- 

v(e* (x 0 ) - e(x 0 )) = e(e (x 0 ) - e{x 0 ))" 



= H ^,^j E (Mx, )e(xj )) - A i E{e{x i )e{x 0 )) + E{e{x 0 )) 2 

i=i j =i 1=1 


usando: 

i. £[f(x 0 )] = 0 

ii. e{ £ (x, )e ( Xj )) = COV(e(x i ), £ {x j )) = C, 

iii. E{ £ {x 0 )) 2 =(7 2 

V( £ *{x 0 )- £ {x 0 ))= -2^^ +a 2 (7) 


derivando respecto a Xi se tiene: 


dv(e' tv,,) - f(.r„)) 3 


+ 2A I £Z J C lJ + - 2A,C W - 2 j>,C, 0 + a‘ 


i—2 7=2 


= 2A.C,, + 2- 2C 10 

7=2 

= 2%A,C v -2 C 10 
1=1 

igualando a cero 

Z2,C„ = C,„ . 

1=1 


En general para cualquier i, i= 1, 2, n, se obtiene: 



EVr c 


7=1 


iO 


Con las /I ecuaciones resultantes se construye el siguiente sistema de ecuaciones: 


( 8 ) 
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f c 
'-'11 

c 

^12 

.. c \ 

'"'l n 

fV 


(c \ 

'-'10 

C 2 i 

^22 

•• c 2 „ 

K 

= 

^20 

S'nX 

C„ 2 • 

Cnn y 

AJ 


l^nOy 


• Varianza de Prediction Kriging Simple. 

Se tiene de (7) que: 

v(e(x„)-e(x„))= - 2 Z^,C,o + o- ! 

i=l j=l i =1 

<r\ = ZAZ^C s -2Z^,C„, + cr 2 

i=l j= 1 i=l 

reemplazando (8) en (7) 

<V = XlCn-lZWn + a* 

i =1 i =1 

a k = ^ 2 -Z^c, 

i=l 


4.4.2. Kriging en Bloques. 

En los dos metodos kriging hasta ahora descritos el objetivo ha estado centrado en la 
prediccion puntual. A menudo, sin embargo, se requiere estimar un bloque, o mas 
precisamente, estimar el valor promedio de la variable dentro de un area local. 



El valor promedio dentro del bloque es estimado por : 


z(A) = £a,.zU i .) 

i =1 
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Del sistema de ecuaciones para el kriging ordinario se tiene: 


'C n • 

• c u r 

'V 


(c ^ 

L 10 

C„i 

l 1 • 

C n „ 1 

• 1 oj 

K 


C„o 

l 1 J 


Consecuentemente el vector del lado derecho de la igualdad en el sistema de arriba debe 
modificarse para incluir las covarianzas respecto al bloque. La covarianza de un punto al 
bloque corresponde a la covarianza promedio entre el punto muestreado i y todos los puntos 
dentro del bloque (en la practica un enmallado regular de puntos dentro del bloque es usado 
como se muestra en la figura de la pagina anterior). El sistema de ecuaciones del kriging en 
bloques esta dado por: 


( r 
^11 

• • C ln 

r 

'V 


(c \ 

c n] • 

C nn 

i 




1 1 ■ 

. . 1 

0, 



l 1 j 


donde el vector de covarianzas al lado derecho de la igualdad en el sistema anterior es 
contiene las covarianzas entre las variables Z(x 1 ), Z(x 2 ), • • ■. Z (x n ) y el bloque A donde se 
quiere hacer la estimacion. 


C„ =Alc„. 

H jljeA 


La varianza del error de prediccion del kriging en bloques esta dada por: 

f n \ ^ 

X X igual a la covarianza entre 


®kB ~ C a A 


+ a 

V .-=1 J 


con Caa = 


w 


Hie.A j! jeA 


pares de puntos dentro del bloque. 

Isaaks y Srivastava (1989) muestran a traves de ejemplos que el kriging en bloques 
coincide con el promedio de predicciones hechas por kriging ordinario sobre cada uno de 
los puntos del enmallado dentro del bloque. Asf mismo indican que en la practica es 
suficiente con un enmallado cuadrado (6x6) para obtener estimaciones estables en los 
bloques. 
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4.4.3. Kriging Universal. 

En los supuestos hechos hasta ahora respecto a los metodos kriging se ha asumido que 
la variable regionalizada es estacionaria (al menos se cumple con la hipotesis intrfnseca). 
En muchos casos, la variable no satisface estas condiciones y se caracteriza por exhibir una 
tendencia. Por ejemplo en hidrologla los niveles piezometricos 4 de una acuffero pueden 
mostrar una pendiente global en la direccion del flujo (Samper y Carrera, 1990). Para tratar 
este tipo de variables es frecuente descomponer la variable Z(x) como la suma de la 
tendencia, tratada como una funcion determimstica, mas una componente estocastica 
estacionaria de media cero. Asuma que: 

Z(x) = m(x) + £(x) 

con E{e(x)) = 0, V{e{; x)) = a 2 y por consiguiente E(z(x)) = m(x). 

La tendencia puede expresarse mediante: 

p 

m[x)=Y J a ,fi{x) 

1 =1 

donde las funciones f,{x) son conocidas y p es el numero de terminos empleados para 
ajustar m(x). 

El predictor kriging universal se define como: 

z *(*o ) = i^Z(xi) 

1=1 


este sera insesgado si: 
e{z(x 0 ))= m(x 0 ) 

^fZ^' Z ( x <)] = m ( x o) 


V «=i J 

^ n A 

»•) = m ( X 0 ) 

V 1=1 J 

Z^/Z/U',) = Yj a if,{x 0 ) 

i=i v 1=1 y /=i 

Z a i =Tj a ifi( x o) 

1=1 v 1=1 J 1=1 


Z z ,/,(a)=Zz'/U,) 


4 


Piezometro: Instrumento utilizado para medir coeficientes de compresibilidad de solidos, lfquidos y gases 
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La obtencion de los pesos en el kriging universal, analogo a los otros metodos 
kriging, se hace de tal forma que la varianza del error de prediccion sea minima. 


v(z-(x 0 ) - Z(x „)) = e(z- (x 0 ) - Z(x 0 )f 


A 


A 2 

) 

J 


= M*, ) “ £ ( X i )) - M*o) - £ ( x o )) 

f V f \ 

X ^M x i ) - m (- x o ) +| X ^ £ ( x i ) - £ ( x o ) 


“|2 


= E 


\i =1 
f 


X^)-^o) 

V «=i J 


J V »=.i 

A 2 " 


J 


= X X 2,2 jE( £ ( x i )e(xj )) - 2X A i E{s{x i )e(x 0 )) + E(e(x 0 )) 2 


Usando 

C„ =COV(e(x,),e(x,)) 

a 2 = E(e(x n )f 
se tiene 

vfz’Uo)- z(xj= X tw,c, ~ 2 i^ c , 

1=1 j=l 1=1 

Luego incluyendo la restriction dada por la condition de insesgamiento, se debe 
minimizar: 

^ 'XV.W-zW 


+ a 2 . 


°lu =XZW® - + o- 2 + 2>, 

1=1 7=1 i =1 /=1 

o en terminos de la funcion de semivarianza 


/=i 


+ 2 X4r io + Xa 

|'=1 7=1 


i=l 


/=1 


X^//( x /) _ Z( x o) 


derivando la expresion anterior respecto a A n A 2 , ■■■, A n , ju x , ju 2 ,---, jU p e igualando a cero 
las correspondientes derivadas se obtienen las siguientes ecuaciones: 


n r 

+ TjVifi (*,) = Lo i = A 2,« 

r=i /=i 

X 2,// (xj ) = /, (*„ ) j = 1, 2, p 


7=1 
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en terminos matriciales 


>11 

7n 

••• Yu 

fn 


f v 


r Y w ' 

7n 

722 

••• 7 ln 

fn 

C4 



720 

YnX 

7n2 

••• Yn2 

fu 

... f 

«z pn 


= 

Yn 0 

fu 

fn 

••• fu 

0 

... o 

A, 


fw 

Jpl 

f P 2 

••• fpn 

0 

... o y 

^nJ 


>0/ 


donde f tj = f l (vjes la 1-esima funcion en el punto j-esimo. 

La varianza de prediccion del kriging universal esta dada por (Samper y Carrera, 1990): 

= Xf-7 l0 + XufiUo)- 

i=i 1=1 

Notese que si p = 1 y f,{x) = 1, el sistema de ecuaciones del kriging universal y la varianza 

de prediccion coinciden con las del kriging ordinario. En este orden de ideas puede decirse 
que el kriging ordinario es un caso particular del kriging universal. 

4.4.4. Kriging Residual. 

La tecnica kriging residuales empleada bajo las mismas circunstancias del kriging 
universal, es decir en aquellos casos en que la variable regionalizada no es estacionaria 
debido a la presencia de tendencia espacial en el valor promedio de la variable. La hipotesis 
central del kriging residual consiste en suponer conocida la tendencia m(x). A partir de ella 
se calculan los residuos con base en los cuales se aplica kriging ordinario. La estimacion de 
la tendencia es generalmente llevada a cabo por medio de mfnimos cuadrados. La 
prediccion en un sitio no muestreado es igual a la tendencia estimada mas la prediccion del 
error, es decir: 

Z (x 0 ) = m(x 0 ) + e (v 0 ) 

e *(*o)=Xv(*i) 

i=1 

los pesos o ponderaciones son estimados por kriging ordinario como se muestra en la 
seccion 4.2. La varianza de prediccion de la variable de interes coincide con la varianza de 
prediccion de los errores. En la figura 15 se muestra un esquema con el procedimiento 
kriging residual en el caso de una tendencia lineal. 
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4.4.5. Kriging Indicador 


Suponga que se tiene una variable regionalizada |z(x): xe D cz R d ] estacionaria. Se 
define la siguiente transformacion: 

J( 7U Ji SiZ( Xi )< Zl 
' [0 Otro caso 

Algunas propiedades son las siguientes: 

i. Pr(/(x, z,) = l) = Pr(z(x) <z,) = F{z ,) 

ii. E(l(x, z, )) = lPr(/(x, z, ) = l) + 0Pr(/(x, z,) = 0) 

= lPr(/(x, z,) = l)= F(z,) 

El predictor kriging indicador es igual a: 

n 

Z l ) = Y J ^Xz l )l(x i , Zi) 

i =1 

es decir que la prediccion de la funcion indicadora en el sitio xo es igual a una combinacion 
lineal de las n funciones indicadoras evaluadas en los sitios de medicion. Samper y Carrera 
(1990) muestran que el kriging indicador es un estimador de la probabilidad acumulada 
hasta el lfmite z definido en la funcion indicadora. El predictor kriging indicador (dado que 
predice probabilidades acumuladas) tiene las siguientes restricciones: 

i. 0< I*{x, z,)< 1 

ii. I*{x, Z[ )< / (x, z'i) si z l < z, 

Una condicion suficiente para que estas restricciones se cumplan es que 
A, i (z l ) = A,, con 0 < A, < 1, Vi, Vz,. 

Sin embargo en la practica las ponderaciones se estiman de tal forma que el predictor sea 
insesgado de varianza minima. 

Para la condicion de insesgamiento: 

E(l* (*o> z ,)) = E(l(x 0 ,z,))= F(z ,) 

Y J ^ l (z,)E{l{x l ,z l ))=F(z,) 

i =1 

J/l i{z,)F{z,)= F{z,) => YjA,{z,) = 1 

i =1 i=l 

Despues de llevar a cabo el proceso de derivacion sobre la expresion de la varianza del 
error de prediccion (obtenida de forma analoga a como se hizo en kriging ordinario), se 
obtiene el siguiente sistema de ecuaciones: 
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n 

Z 4- (*/ K + P = r,o i =1 ’ 2, ...,n 
1=1 
n 

i=l 

donde 7, 7 = /(/;), la funcion de semivarianza evaluad para la distancia entre los sitios i, j. 

La varianza del error de prediccion se encuentra de la misma forma a como se menciono en 
la seccion 4.2. 


4.4.6. Kriging Log-Normal y Multi-Gaussiano 


En estos dos procedimientos se hacen transformaciones de la variable regionalizada 
con el proposito de normalizar en cada sitio de la region de estudio. 

El primero de estos consiste en aplicar kriging ordinario a la transformacion logaritmica de 
los datos. Veamos: 

Sea {z(x):xeZ)} una variable regionalizada log-normal. Es decir que 
F(x) = Log(z(x)) tiene distribucion normal. Algunas veces se requiere adicionar una 
constante positiva de tal forma que Y(x) este definida. 

El predictor kriging log-normal es: 

rK)=jv(i,)- 

1=1 

Los pesos se obtienen de manera analoga al kriging ordinario. El semivariograma 
usado es el de los valores transformados. La complicacion practica puede darse al hacer la 
retransformacion a la escala original, puesto que Z (x 0 ) = cxpfk" (x 0 )) es un predictor 


sesgado. 

Se puede demostrar que un predictor insesgado es (Cressie, 1993): 


Z'(x 0 )=exp rW^-4 dende <7;„es la varianza de prediccion obienida en el 


sitio xo por medio de kriging ordinario sobre los valores transformados y // es cl 
multiplicador de Lagrange empleado para la condicion de insesgamiento sobre la escala de 
valores transformados. 

Respecto al kriging multi-gaussiano, suponga que se tiene una variable 
regionalizada {z(r):ie D} estacionaria. Este procedimiento consiste en hacer una 
transformacion de Z(x) tal que los valores transformados sigan una distribucion normal 
estandar. En ese sentido es una generalizacion del kriging log-normal. 

Los pasos del metodo kroging multi-gaussiano son los siguientes: 

i. Se encuentra la funcion de probabilidad acumulada empfrica F n (z(x)). 

ii. Se calculan con base en F n (z(x)) los “scores” normales estandar (Fig. 16), es decir 
los valores de una distribucion de probabilidad normal estandar para los cuales la 
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probabilidad acumulada corresponde a F n (z(x)). En otras palabras se encuentra 

i/W=*-'(^(zW)). 



z ( x ) U(x) 

Figura 16. Representacion de la transformacion a scores normales 


Hi. Se realiza kriging simple sobre los valores transformados. 


4.5. Aplicacion : Estudio de la distribucion espacial de variables fisicoquimicas y 
bioldgicas medidas en el estuario Cienaga Grande de Santa Marta 

Con base en la informacidn descrita en la aplicacion de la seccion 3.3. y empleando 
los semivariogramas ajustados a las variables alii mencionadas, se generaron mapas de 
isolmeas de cada una de ellas (figuras 17 a 22) y se realizo interpretacion de los mismos en 
un contexto ecologico. Como apoyo en la descripcion se emplean medidas descriptivas de 
las variables (tabla 6). 


Tabla 6. Medidas descriptivas de variables fisicoquimicas y biologicas medidas durante un muestreo 
realizado en marzo de 1997 en la Cienaga Grande de Santa Marta. 


Variable 

Nivel de la 
columna 

Promedio 

Minimo 

Maximo 

Coeficiente de 
Variacion (%) 

Salinidad 

Superficie 

17.6 

13.02 

34.9 

16.1 


Fondo 

18.04 

13.94 

33.9 

15.5 

Oxigeno 

Superficie 

8.80 

3.03 

16.2 

36.9 


Fondo 

5.68 

2.64 

13.4 

36.8 

Solidos en suspension 

Superficie 

218.2 

103 

318 

18.8 


Fondo 

215.35 

86 

310 

19.6 

Nitritos 

Superficie 

0.43 

0.01 

1.61 

70.8 


Fondo 

0.42 

0.01 

2.39 

81.7 

Clorofila a 

Superficie 

132.44 

2.91 

198.35 

23.8 


Fondo 

136.19 

2.91 

194.75 

26.4 

Profundidad 


1.47 

0.25 

2.50 

24.1 

Transparencia 


0.27 

0.20 

0.35 

10.8 
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4.5.1. Mapas de Distribution 
Salinidad. 

Los valores medidos oscilaron entre 13.02 y 34.9 en el nivel superficial de la columna 
de agua y entre 13.94 y 33.9 en el fondo de la misma, con valores promedios de 17.6 y 
18.04, respectivamente (tabla 6). De los mapas de isolmeas de distribucion de dicha 
variable (Fig. 17), es posible concluir que existe gran homogeneidad en todo el cuerpo de 
agua, con excepcion de la zona nororiental, puesto que los valores maximos y mfnimos 
predichos varfan solo alrededor de 14 y 19 unidades. 

En los dos mapas (superficie y fondo) de la figura 17, se evidencia la influencia que 
tienen las entradas de agua dulce y marina sobre la magnitud de la variable dentro del 
sistema. Las salinidades maximas se encuentran en la zona nororiental (zona estuarina) 
donde hay entrada de agua marina a traves del sitio denominado Boca de la Barra (Fig. 1), 
encontrandose alii valores superiores a 30 unidades. Hacia la zona centro del cuerpo de 
agua se presentan las menores magnitudes de la variable, valores entre 15 y 16 unidades, lo 
que parece ser consecuencia del aporte de agua dulce que se da en la desembocadura de uno 
de los tres rfos (Rio Sevilla) que baja de la Sierra Nevada de Santa Marta (SNSM). Asi 
mismo se puede observar que en el sector occidental del sistema se presentan valores 
intermedios a los de las zonas antes mencionadas (alrededor de 19 unidades). Lo anterior 
puede deberse al efecto de intercambio de aguas, por medio de los canales Grande y Clarin, 
con el ecosistema Complejo Pajarales (Fig. 1), en donde se da un proceso de lavado de 
suelos hipersalinos en epoca de lluvias o cuando hay inundaciones. Dada la similitud en 
magnitud y forma de distribucion que se observa en los mapas de superficie y fondo (Fig. 
17), se podria pensar que para la epoca seca del ano, no existe estratificacion de la columna 
de agua respecto a esta variable. 

Los valores de salinidad observados y predichos a traves de la tecnica kriging, resultan 
bajos respecto a los registrados para esta misma epoca en otros estudios (Giraldo et al., 
1995). Lo anterior podria deberse a un posible aumento en los caudales de los rios que 
desembocan en la CGSM, durante los meses de lluvia precedentes al muestreo, como 
consecuencia del efecto del fenomeno del nino en la region a finales del ano 1996. No 
obstante lo anterior, puede pensarse, dada la gran homogeneidad en la distribucion, que 
para la fecha del muestreo no se estaban presentando aportes considerables de agua dulce, 
por parte de los rios que desembocan en la CGSM, lo que significa un periodo de relativa 
calma para los organismos que dependen de la salinidad para sus funciones y distribucion 
(Reidy Wood, 1976). 

Oxigeno Disuelto. 

Los valores medidos de esta variable oscilaron entre 3.03 (mg/1) y 16.2 (mg/1) en la 
superficie de la columna de agua y entre 2.09 (mg/1) y 13.4 (mg/1) en el fondo de la 
misma, con valores promedios de 8.8 (mg/1) y 5.68 (mg/1), respectivamente (tabla 6). 

Las correspondientes isolmeas (Fig. 18), indican que en el fondo de la columna de agua se 
presenta mayor homogeneidad en la distribucion, puesto que los valores predichos varfan 
entre 4.5 mg/1 y 6.5 mg/1, con excepcion de una pequena zona en el sector nororiental frente 
a la desembocadura del rfo Sevilla (valores entre 6.5 y 9.5 mg/1), mientras que en superficie 
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existe considerable diferencia entre los valores ajustados en el centro del sistema (entre 9 
mg/1 y 13 mg/1) y los estimados para la zona sur y noroccidental del mismo (magnitudes 
alrededor de 4 mg/1). Lo anterior sugiere la ocurrencia de procesos de estratificacion en el 
sistema hacia la zona central del espejo de agua, donde la productividad se concentra 
aportando grandes volumenes de oxrgeno al agua durante el dra (Reid y Wood, 1976; 
Welch, 1992; Mancera y Vidal, 1994). Las isolmeas, para ambos niveles de la columna de 
agua, muestran que hacia las fronteras del sistema los valores del gas disminuyen. Este 
comportamiento podrfa ser explicado al considerarse que en estas zonas existe intercambio 
de flujos entre el sistema y otros cuerpos de agua, ademas de aportes de hojarasca y 
material organico, provenientes del manglar. 

Solidos en suspension. 

Los valores para la variable, presentan algunas diferencias entre los dos pianos de 
muestreo. Los mrnimos y maximos fueron de 103 mg/1 y 318 mg/1 en la superficie de la 
columna de agua y de 86 mg/1 y 310 mg/1 en el fondo de la misma, con promedios de 218.2 
mg/1 y 215.3 mg/1, respectivamente (tabla 6). 

El mapa de distribucion superficial (big. 19), revela la influencia que tienen los aportes de 
agua sobre la magnitud de esta variable en el sistema. Se observa que las mayores 
concentraciones se presentan en las zonas de las desembocaduras de los rros Fundacion y 
Aracataca ademas de la del cano Clarfn (por medio del cual se da el aporte de agua del rfo 
Magdalena) y que las menores magnitudes se dan en el sector de intercambio de agua dulce 
y marina (desde la zona centro y nororiental hacia el sitio denominado Boca de la Barra). 
Una excepcion a este comportamiento se da en la desembocadura del rfo Sevilla y del cano 
Grande en donde las concentraciones de los solidos en suspension son muy similares a las 
observadas en el resto del cuerpo de agua. 

El patron de comportamiento de la variable en el fondo de la columna de agua es muy 
similar al descrito en el parrafo de arriba respecto a los valores superficiales; es decir 
mayores concentraciones hacia las desembocaduras de los rros y canos (zonas sur y nor¬ 
occidental) y menores magnitudes en la zona nororiental. Sin embargo, la diferencia entre 
los valores predichos en estas fronteras y los del resto del sistema (valores entre 220 mg/1 y 
210 mg/1), no resultan significativos como en el caso de la distribucion superficial (valores 
entre 245 mg/1 y 175 mg/1). 

Una posible explicacion a la diferencia en magnitud de los valores de superficie y 
fondo es que los flujos de agua dulce son menos densos y presentan mayores 
concentraciones de solidos en suspension, por lo cual al ingresar al sistema y encontrarse 
con las aguas salobres del mismo (mas pesadas), tienden a permanecer en la superficie (lo 
cual puede causar estratificacion. (Wheaton, 1977; Welch, 1992 ; Jay et al, 1997; Nixon, 
1997). 

Nitritos. 

El ion nitrito presento valores entre 0.01 p mol/1 y 1.61 pmol/1 para la superficie de la 
columna de agua y entre 0.01 pmol/1 y 2.39 pmol/1 en el fondo de la misma. Los valores 
promedios fueron de 0.43 pmol/1 y 0.42 pmol/1, respectivamente (tabla 6). 
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Las mediciones superiores a 1 pmol/1 se dieron en dos estaciones de muestreo, ubicadas en 
la zona norte del sistema (cuadrfculas 14 y 15, Fig. 1). Dado que lo anterior no fue el patron 
generalizado, las isolmeas de distribucion en superficie y fondo (Fig. 20) presentan solo 
valores alrededor de los promedios arriba mencionados. Este resultado es esperable, puesto 
que los nitritos generalmente se dan en bajas concentraciones (Mancera, 1990; Hernandez y 
Gocke, 1990). Day et al. (1989), indican que esto puede ser debido al consumo continuo de 
las comunidades fitoplanctonicas y a la precipitacion en los sedimentos como consecuencia 
de lo cambios en las condiciones del agua estuarina. 

En ambos casos (superficie y fondo) los valores interpolados (entre 0.2 pmol/1 y 0.7 
pmol/1) para esta variable, revelan la presencia de un gradiente sur-norte, dandose las 
mayores concentraciones en el sector mas estuarino. Lo anterior podrfa ser consecuencia de 
aportes de materia organica por parte de las poblaciones cercanas a esta zona (Welch, 
1992). Los mapas de distribucion espacial no revelan estratificacion de la columna de agua 
para esta variable, dada la similitud en los valores predichos en superficie y fondo (Fig. 
20 ). 

Clorofila “a”. 

Los valores medidos de clorofila “a” oscilaron entre 2.91 pg/1 y 198.35 pg/1 en la 
superficie de la columna de agua y entre 2.91 pg/1 y 194.75 pg/1 en el fondo de la misma. 
Los valores promedios fueron 132.44 pg/1 y 136 pg/1, respectivamente. Los bajos 
coeficientes de variacion (menores del 30%), en ambos casos (superficie y fondo), indican 
relativa homogeneidad en las mediciones de esta variable (tabla 6). Las isolmeas de 
distribucion calculadas con los datos predichos (Fig. 21) presentan algunos aspectos 
comunes. En ambos mapas (superficie y fondo) se observa que los valores maximos 
(alrededor de 160 pg/1) se dan en el sur del sistema hacia la desembocadura del rfo 
Fundacion y las menores concentraciones se presentan en el sector mas nororiental (valores 
menores de 50 pg/1). La diferencia radica en el comportamiento en la zona centro del 
espejo de agua. Mientras que en la superficie se presenta alta variabilidad (valores entre 90 
y 160 pg/1), en el fondo de la columna de agua los valores son muy homogeneos (entre 130 
pg/1 y 150 pg/1) y se ajustan claramente a una tendencia creciente en sentido sur - 
nororiente. El comportamiento distribucional representado en los mapas de isolmeas puede 
estar de acuerdo con las condiciones climaticas de la epoca. Bula-Meyer (1989) y Sanchez 
(1996), indican que en la epoca mas seca del ano (epoca en la que se realizo el muestreo) 
predominan los vientos Alisios y que la circulacion de las masas de agua en el sistema 
obedece a la fuerza del viento y a los cambios de marea en el Mar Caribe. Esto hace pensar 
que los flujos de agua son mas lentos en la zona centro del sistema y por consiguiente, al no 
presentarse un recambio de agua muy fuerte, se favorece el desarrollo de las comunidades 
fitoplanctonicas, puesto que estas consumen los nutrientes que se liberan desde el 
sedimento por accion de los vientos (Welch, 1992). 

Giraldo (1996), encontro un comportamiento similar en la distribucion de esta variable 
con datos promedios de la epoca, calculados con informacion de varios anos anteriores a 
1995. Sin embargo en la zona de la desembocadura de los rfos, especfficamente en la del 
Fundacion, los valores reportados por dicho autor resultaron considerablemente mas bajos a 
los encontrados en el presente estudio. Lo anterior puede estar indicando, como se 
menciono en la interpretacion de los resultados obtenidos con la salinidad, una disminucion 
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de la entrada de agua dulce en esta zona, causandose asf el mismo efecto de baja 
circulacion y alta asimilacion de nutrientes, comentado en el parrafo de arriba para la zona 
centra del sistema. 

Profundidad (Batimetria) y Transparencia (Secchi). 

Las medidas resumenes de la variable profundidad (tabla 6) y el mapa de distribucion 
calculado con los datos predichos (Fig. 22) confirman lo reportado en estudios anteriores 
(Wiedemann, 1973), en los cuales se afirma que la CGSM es un sistema somero, con una 
profundidad promedio alrededor de 1.5 m. El mencionado mapa revela la presencia de 
gradientes positivos en sentido oriente occidente y sur nor-occidente, respectivamente. Este 
comportamiento puede estar relacionado con los procesos de sedimentacion que se han 
venido presentando en los ultimos anos cerca al sitio Boca de la Barra (cuadrfcula 1, Fig. 
1) como consecuencia de la disminucion de los flujos de agua que entran al sistema a traves 
de rfos que bajan de la SNSM y de los canales que comunican con el no Magdalena. 

De otro lado respecto a la variable transparencia, medida por medio de la profundidad 
del disco de Secchi (Reid y Wood, 1976) se puede concluir, de acuerdo con los valores 
encontrados (tabla 6) y predichos (Fig. 22), que esta es una caracterfstica muy homogenea 
en el sistema de estudio. Fas isolmeas de distribucion indican que en general los valores 
esperados en este ecosistema no son superiores a 30 cm. Este bajo nivel de transparencia 
esta de acuerdo con el resultado encontrado para la variable clorofila “a” (altas 
concentraciones en gran parte del sistema, como consecuencia de un aumento en las 
entradas de nutrientes), puesto que como lo muestra Welch (1992), la relacion entre estas 
dos variables es de tipo inverso. Segun resultados reportados por este autor se espera que 
para niveles de visibilidad del disco secchi, inferiores a 1 m se den concentraciones de 
clorofila “a” superiores a 80 pg/1. 
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Figura 17. Distribucion espacial de la salinidad del agua en la Cienaga Grande de Santa Marta durante una jornada de muestreo realizada en marzo de 1997. El mapa de 
la izquierda corresponde a los valores en la superficie de la columna de agua y el de la derecha a los niveles en el fondo de la misma. 
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Figura 18. Distribucion espacial del oxfgeno disuelto (mg/1) en la Cienaga Grande de Santa Marta durante una jornada de muestreo realizada en marzo de 1997. El mapa 
de la izquierda corresponde a los valores en la superficie de la columna de agua y el de la derecha a los niveles en el fondo de la misma. 
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Figura 19. Distribution espacial de los solidos en suspension (mg/1) en la Cienaga Grande de Santa Marta durante una jornada de muestreo realizada en marzo de 1997. 
El mapa de la izquierda corresponde a los valores en la superficie de la columna de agua y el de la derecha a los niveles en el fondo de la misma. 
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Figura 20. Distribution espacial de nitritos (|imol/l) en la Cienaga Grande de Santa Marta durante una jornada de muestreo realizada en marzo de 1997. El mapa de la 
izquierda corresponde a los valores en la superficie de la columna de agua y el de la derecha a los niveles en el fondo de la misma. 
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Figura 21. Distribucion espacial de clorofila a ((J.g/1) en la Cienaga Grande de Santa Marta durante una jornada de muestreo realizada en marzo de 1997. El mapa de la 
izquierda corresponde a los valores en la superficie de la columna de agua y el de la derecha a los niveles en el fondo de la misma. 
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Figura 22. Distribucion espacial de la profundidad (m) (izquierda) y transparencia (m) (derecha) en la Cienaga Grande de Santa Marta durante una jornada de muestreo 
realizada en marzo de 1997. 
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Capitulo Cinco 


Temas Especiales 

En este capitulo se presentan algunos temas que no aparecen explfcitamente en libros 
clasicos de geoestadfstica, tales como el diseno de redes muestrales o el analisis de 
componentes principales sobre variables regionalizadas. Asf mismo se revisa la teorfa del 
analisis cokriging y de simulacion de fenomenos espaciales, bajo el supuesto de normalidad 
multivariada. 

5.1. Cokriging Ordinario 

Si se tienen dos variables regionalizadas Z v i(x) y Z V 2 (x) tomadas en sitios de la region de 
estudio, no necesariamente iguales, entonces el semivariograma cruzado entre ellas, se 
estima por: 

1 nh 

n,v 2 w=- — 2 ] i z v, & + h )~ z v, (*w V2 y+ h )~ z V2 (*)} (9) 

Zn h 

Donde n/, es el numero de parejas de datos que se encuentran a una distancia h (Bogaert et 
al, 1995). 

• Modelo Lineal de Corregionalizacion (MLC) 

El MLC asume que todos los semivariogramas simples y cruzados pueden expresarse 
como una suma de modelos basicos (exponencial, esferico, gaussiano, etc.) identicos. Para el 
caso de dos variables: 


Tv x ( h)=a 0 y 0 ( h)+...+a m y m (h) 
Tv 2 (h)^ Q r Q (h)+... +J 8 in r m (h) 
Yv x v 2 ( K)=5qYq (h)+...+S m y m (/?) 


( 10 ) 


donde y v /(h) y yMh) son los semivariogramas simples, y v i v 2 (h) es el semivariograma cruzado. 
yo(h), y(h), . . yjh) son los modelos basicos de semivariograma y a, P y 8 son constantes. 


Matricialmente: 
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Y s (h) 


( 11 ) 


IJh) se le conoce como matriz de corregionalizacion. 
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• Predictor Cokriging 

El metodo de prediccion espacial cokriging consiste en hacer prediccion espacial de una 
variable con base en su informacion y en la de algunas variables auxiliares que este 
correlacionadas espacialmente con ella. El predictor cokriging tiene la siguiente expresion en 
el caso en el que se considera una sola variable auxiliar: 

n \ 

K (Xj) (i2) 

i =1 7=1 

el lado izquierdo de la igualdad en la ecuacion anterior representa la prediccion de la 
variable de interes en el punto x 0 no muestreado. Z V Jx t ) con i=l, 2 , ... , it/, representa la 

variable primaria. As! mismo, Z v J x ] ) con j=l, 2, n 2 , representa la variable auxiliar. a, 

y bj, con i-1, 2 , , m y j=l, 2, .. ., n 2 respectivamente, representan los pesos o 

ponderaciones de las observaciones de las variables primaria y auxiliar y se estiman con base 
en el MLC ajustado a los semivariogramas simples y cruzados. Los pesos a, y bj se estiman 
de manera analoga al proceso descrito para el metodo kriging ordinario, es decir estos seran 
los que minimizan la varianza del error de prediccion sujeta a la restriccion de que el 
predictor sea insesgado. La estimacion de los parametros se obtiene resolviendo el siguiente 
sistema de ecuaciones (Isaaks y Srivastava, 1989): 

' %m - rJni) r,JW - 1 oVaA r ,m' 

yjln) ••• yjt%n) y.JIn) ••• y vlv2 (mn) 1 0 a n yJQn) 

y,JV) ••• yjul) yjll) ••• yjml) 0 1 b, y,JO,l) 

; ; ; ;;;;;; ; (i3) 

y,Jl,tn) ••• yjivn) yjlm) ••• yjmm) 0 1 b m y lv2 (0,m) 

1 ••• 1 0 ■■■ 0 0 0 m 1 

v o ••• o i ••• i o oJUJ l 0 , 

La matriz del lado izquierdo contiene los valores de las funciones de semivarianza y de 
scmivarianza cmzada calculadas para todas las distancias entre las parejas de puntos 
consideradas. Las dos ultimas filas de dicha matriz son las correspondientes a la restriccion 
de insesgamiento del predictor, a, y bj con i = 1 , 2, ..., n y j = 1 , 2, ...., m, son los parametros 
a estimar. pi y p 2 son los multiplicadores de Lagrange empleados para la restriccion de 
insesgamiento y el vector del lado derecho contiene los valores de la funciones de 
semivarianza y semivarianza cmzada evaluados para las distancia entre los sitios de 
muestreo (de ambas variables) y el sitio donde se desea hacer la prediccion. Las dos ultimas 
filas del vector estan asociadas a la condicion de insesgamiento. La correspondiente 
varianza de prediccion del metodo cokriging se calcula como (Bogaert et al, 1995): 

a\ = Co\{Z vI (x 0 ), Z vI (x 0 )) + //, +!>, Cov(Z vl (x ,), Z vl (x 0 )) + £ b } Cov[z v2 {x J ), Z v2 (x 0 )) ( 11 ) 

i=i i=i 
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donde |ii es el multiplicador de Lagrange empleado para la restriction dado por la condition 

( n 

de insesgamiento ^ J a i =l . 

\i—l J 

Cov(Z vi (x k ), Z vi (x ,)) = a 2 vl - Y vlvi {k,l) es la funcion de covarianza espacial de la variable i, 


i=l,2, evaluada para la distancia entre los sitios de muestreo k, l. 

La ventaja del metodo cokriging frente al kriging radica en el hecho de que cuando la 
variable auxiliar esta ampliamente correlacionada con la variable de interes se puede obtener 
un diminution en la varianza de prediccion, no obstante dicha variable tenga menor 
densidad de muestreo. En situaciones en las que la valuable objetivo tiene costos altos de 
muestreo se recomienda la aplicacion de esta metodologia (Bogaert et al., 1995). 


• Kriging Probabilistico 


Es un predictor basado en cokriging que utiliza como valuables predictoras una valuable 
indicadora y una valuable generada a traves de la transformation uniforme. 

Sea Z(xi) la variable observada, i = 1,2, . . ., n, entonces se definen las siguientes 
transfoiunaciones: 








1 SiZ(Xj)<z 
0 Otro caso 


U (x ,) 


R(Z( Xi ) 

n 


para todo i, i = 1,2,..., n. 


con R{Z(xi )) igual al rango (posicion que ocupa dentro de los datos ordenados de menor a 
mayor) la i-esima observacion muestral. La prediccion de probabilidad de exito en el sitios 
de interes esta dada por: 


n n 

1 (*o)= X (x ' ’ z)+ X v i u(x i ) 

i=i (=i 

Los pesos y Vi se estiman mediante el sistema de ecuaciones del metodo cokriging. 

5.2. Analisis de Componentes Principales Regionalizado (ACPR) 

El ACPR se fundamenta en la realizacion de analisis de componentes principales (ACP) 
(apendice, seccion 6.4.2) con base en varias matrices de corregionalizacion (seccion 5.1). El 
caso mas simple de ACPR es cuando se aplica ACP con base en la matriz de 
corregionalizacion a distancia cero (matriz de correlacion tradicional). En este caso la 
tecnica consiste en generar los ejes principales de la forma tradicional (Manly, 1994), 
posteriormente realizar la correspondiente interpretacion de estos en terminos de la 
variabilidad explicada por cada componente respecto a cada valuable original y finalmente 
llevar a cabo un analisis geoestadfstico a traves de la estimacion de la funcion de 
semivarianza y de la aplicacion de algun procedimiento kriging con base en los datos de los 
ejes generados. La interpretacion del mapa de predicciones obtenida sobre los componentes 
permite obtener una vision integral del comportamiento conjunto de las variables 
consideradas dentro del sistema de estudio. 
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El procedimiento se puede realizar con matrices de correlacion obtenidas para diferentes 
distancias (fijadas de antemano) y calculadas a traves de la funcion de semivarianza, 
covarianza cruzada o de correlacion cruzada (si las variables tienen diferentes escalas de 
medida se recomienda emplear la funcion de correlacion cruzada). Lo anterior implica que 


se deben calcular 



funciones de correlacion espacial cruzada, siendo n el numero total de 


valuables involucradas en el estudio. Esto puede ser una limitante computacional del 
metodo, cuando se incremente el numero de variables. En la practica se acostumbra a 
seleccionar grupos de pocas valuables (alrededor de 5, consideradas como las mas 
relevantes) que esten muy relacionadas espacialmente y con base en la informacion de estas 
hacer el analisis para dos o tres matrices de corregionalizacion (incluyendo la de distancia 
cero). 


• Correlacion Intrinseca. 


Se dice que un conjunto de valuables regionalizadas tiene correlacion intrinseca cuando 
la estructura de correlacion de las valuables es independiente de la distancia espacial (puede 
haber correlacion para distancia cero), es decir cuando las funciones de semivarianza 
cruzada, covarianza cruzada o correlacion cruzada, entre las valuables, son constantes en 
funcion de la distancia. La detection de correlacion intrinseca puede hacerse a traves de las 
valuables originales o por medio de los componentes principales generados. En la section 
6.4.2., se menciona que los ejes principales deben ser independientes, luego se espera que 
para cualquier distancia la funcion de semivarianza cruzada entre cualquier para de 
componentes principales este alrededor de cero. En caso contrario habra relation espacial 
entre las variables. Existen dos formas de llevar a cabo el ACPR dependiendo de si las 
variables tienen o no correlacion intrinseca. 


• ACPR en Presencia de Correlacion Intrinseca. 

El algoritmo en este caso se inicia con el calculo de la matriz de corregionalizacion para 
distancia cero (matriz de correlacion clasica) (en las otras distancias no hay correlacion entre 
las valuables). Posteriormente se aplica el ACP clasico se generan ejes factoriales que 
explican, se espera que en un alto porcentaje, la variabilidad contenida en el conjunto total 
de variables (idealmente dos o tres componentes deberfan explicar mas del 90% de la 
varianza total). Con base en la magnitud y le signo de los vectores propios se identifica el 
peso de cada valuable original en los correspondientes ejes. Finalmente se obtiene un mapa 
de distribution espacial, cumpliendo con las etapas basicas del analisis geoestadfstico, que 
permita dar una interpretation simultanea del comportamiento de las variables involucradas 
en el analisis. 


• ACPR sin Correlacion Intrinseca. 

En este caso se debe establecer en primer lugar un modelo lineal de corregionalizacion 
entre las valuables originales. Con base en este se calculan varias matrices de 
corregionalizacion (una para cada distancia h fijada) y con cada una de ellas se realiza un 
ACP clasico. Los resultados obtenidos en cada analisis peiuniten establecer relaciones entre 
las valuables que no son observados en analisis clasicos de correlacion . 
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5.3. Diseno de Redes de Muestreo 

• Selection de Variables 

Cuando se va a iniciar el estudio de un ecosistema natural, deben establecerse aspectos 
referentes a su ubicacion geografica, a las condiciones climaticas, hidricas y geologicas del 
mismo. La revision debe incluir los antecedentes historicos de las variables ecologicas e 
incluso los factores economicos, socio-culturales y demograficos que enmarcan a la region 
de estudio. Todos estos aspectos permiten planear, desde un punto de vista logistico la 
realization del muestreo. 

Podrfa decirse que existen dos tipos de variables que deben tenerse en cuenta cuando se 
realiza un estudio ambiental. Aquellas que estan directamente relacionadas con el fenomenos 
ecologico de estudio (contamination, calidad del agua, manglar, bentos, plancton, etc.) y 
aquellas que permiten decidir respecto a aspectos relacionados con las frecuencia del 
muestreo (valuables climaticas) y la ubicacion estrategica de sitios de muestreo (por ejemplo 
variables hidricas como direction y velocidad de las masas de agua, fuentes de entrada y 
salida, etc) 

Asi mismo la selection de las valuables es consecuencia del tipo de estudio que se quiere 
llevar a cabo. Estos pueden ser de linea base, monitoreo y evaluation. En el primero de estos 
no hay antecedentes historicos respecto al fenomeno de interes, se asume total 
desconocimiento respecto a la relation, comportamiento y distribution de las valuables en 
el ecosistema y por consiguiente se debe evaluar un mimero grande de variables, con amplia 
intensidad de muestreo en la que se cubra toda la region de estudio, de forma tal que se 
pueda caracterizar de manera general el sistema. Si existe conocimiento de la region de 
estudio, hay estudios preliminares que posibilitan el planteamiento de estructuras de 
correlation espacial y temporal de las variables y se quieren establecer los cambios que se 
estan dando en el ecosistema, por ejemplo por actividades antropicas, es entonces un estudio 
de monitoreo. En este caso debe establecerse con base en la information disponible tanto la 
frecuencia como la ubicacion de los puntos de muestreo. Por ultimo cuando hay 
conocimiento del ecosistema respecto al fenomeno de interes y se quieren observar posibles 
variaciones muy puntuales respecto al patron temporal o espacial tradicionalmente 
observado, el estudio se denomina de evaluation. En este ultimo caso el objetivo puede ser 
el de conservar o mitigar posibles danos mas que el de hacer diagnostico como en el caso del 
monitoreo. 

• Selection de la Red Optima de Muestreo 

Como en cualquier procedimiento estadistico en el que se hace inferencia, en 
geoestadistica cuando se hace prediction en sitios o puntos de la region de estudio no 
observados, a traves de cualquiera de las tecnicas kriging, es necesario evaluar la precision 
de tal prediction. Lo anterior se realiza, como se establecio en el capitulo 4 y en la section 
5.1, por medio del calculo de la varianza del error de prediction. De la section 4.2, para el 
caso del kriging ordinario, la varianza de prediction se calcula por: 


66 



= XAr<» + m 

i =1 


Es evidente, de la ecuacion anterior, que la varianza de prediccion no es constante como 
en el caso clasico y que ademas no depende de los valores medidos de la valuable sino de su 
estructura de correlacion, evaluada a traves de la funcion de semivarianza. 

McBratney et al (1981) muestran que, para cualquier densidad muestral, la distancia 
maxima entre un punto de observacion y un punto a interpolar es minima cuando la 
configuration de los puntos es hecha en un enmallado triangular, por lo cual bajo esta 
distribution de puntos se obtendran las menores varianzas de prediccion. Sin embargo este 
mismo autor y Warrick et al (1986) indican que por razones logrsticas referentes a la 
ubicacion de los sitios en el campo y minimization de recorridos, el enmallado cuadrado es 
preferible. 

De acuerdo con lo anterior el problema del diseno muestral se limita a establecer para 
varias redes de muestreo, de diferentes densidades, con enmallado triangular equilatero o 
cuadrado, la relation entre las varianzas de prediccion maximas (las obtenidas en el centra 
de cada triangulo o cuadrado) y los costos asociados a ellas. De esta forma se deduce el 
costo necesario para alcanzar cierto grado de certidumbre o, a la inversa, la varianza de 
prediccion si se prefija el costo. 

5.4. Simulacion 


A continuation se describe el metodo de simulacion de variables regionalizadas con 
densidad conjunta normal multivariada (gaussiana) (Cressie, 1993) que se fundamenta en la 
tecnica de descomposicion de Cholesky (Anderson, 1984). 


Suponga que se desea simular el vector z(x) 


'zur 

Z(x 2 ) 


correspondiente a n variables 


Uk)J 


aleatorias en n sitios de muestreo de interes x h x 2 , ..., x„. Asuma que el proceso estocastico 
tiene vector de medias y matriz de varianzas y covarianzas 
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V M^X n )y 


S'nl C n2 ’ 

d’nn j 



La matriz de varianzas y covarianzas es descompuesta por el metodo de Cholesky. 

£ = Lll , con L una matriz triangular - inferior. Entonces el vector simulado se define como: 
z(x) = jl{x) + Le , donde e ~ /V„(o, /). 

Usando teoremas referentes a la distribucion de combinaciones lineales de vectores con 
distribucion normal multivariada (Anderson, 1984) se comprueba que el vector simulado 
tiene vector de medias p.{x) y matriz de covarianzas £. 
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5.5. Aplicaciones 

5.5.1. Comparacion de los metodos Kriging y Cokriging con base en resultados de analisis 
espaciales de Variables Fisicoquimicas y Biologicas Medidas en el Estuario Cienaga 
Grande de Santa Marta. 

Se emplea la informacion de las variables profundidad (m), secchi (cm), salinidad, 
oxfgeno (mg/1) y clorofila "a" (pg/1) descrita en la aplicacion 3.3. para mostrar las ventajas 
del metodo cokriging respecto al kriging. 

Se menciono en la seccion 4.1. que la tecnica cokriging es preferible al kriging cuando 
hay informacion de variables auxiliares relacionadas espacialmente con la variable de interes 
y esta ultima tiene altos costos de muestreo. En este caso se usan como variables auxiliares 
profundidad, secchi, salinidad y oxfgeno y como variable objetivo clorofila "a". Las cuatro 
primeras tienen bajos costos de muestreo (son medidas in-situ ) mientras que los altos costos 
de los insumos de laboratorio necesarios para la obtencion de las medidas de clorofila "a" 
pueden ser limitantes del uso de una red de muestreo densa para dicha variable. Aunque no 
se presentan en el trabajo los semivariogramas cruzados entre las variables consideradas 
como auxiliares y clorofila "a", estos muestran la presencia de fuertes estructuras de 
dependencia espacial entre ellas. 

Debido a que la aplicacion del metodo cokriging resulta dispendiosa cuando se 
consideran dos o mas valuables auxiliares (es complejo el ajuste del modelo lineal de 
corregionalizacion), se decidio tomar la informacion de las variables auxiliares de forma 
"condensada" a traves del indicador IGC,(4) (Giraldo, 2002; apendice, seccion 6.1), 
calculado con base en la informacion dicotomizada de dichas variables.. El criterio para 
dicotomizar cada variable fue el de comparacion de cada valor observado con su 
correspondiente mediana. Se asigno el valor 1 para valores mayores o iguales que la mediana 
y 0 en caso contrario. Se aplico el metodo cokriging para realizar predicciones de la valuable 
clorofila "a" con base en sus observaciones e informacion auxiliar del indicador IGC,(4). 
Para detectar la eficiencia del metodo se redujo a aproximadamente la mitad (54 datos) la 
informacion original de clorofila "a" (en adelante se denomina a esta como CLORO 2) y se 
dejo la informacion completa (114 datos) del IGC,(4).E1 ajuste del modelo lineal de 
corregionalizacion y el calculo de las vaiianzas de prediccion fue realizado en el software 
The Spatial Interpolation of Agroclimatic Data (Bogaert et al., 1995). 

• Resultados y Discusion 

En primera instancia se calcularon semivariogramas simples de las variables CLORO 2 
e IGC/(4) y los semivariogramas cruzados entre estas tres variables. Solo se muestra el 
semivariograma experimental cruzado (Fig. 23). Este indica que las dos variables 
consideradas presentan correlacion espacial inversa, es decir que valores altos de 
productividad biologica (alta clorofila "a") estan asociados con valores bajos del indicador 
IGC,(4) en zonas circundantes (incluso mayores a 10 km). Los valores bajos del indicador 
IGC,(4) estan asociados a magnitudes por debajo de la mediana en las valuables profundidad, 
secchi, salinidad y oxfgeno disuelto (ver tabla de interpretacion del IGC,(4) en el apendice). 
De lo anterior se concluye que zonas con alta biomasa fitoplanctonica estan asociadas a baja 
profundidad, alta turbidez y a masas de agua con baja salinidad y bajo nivel de oxfgeno (tal 
vez como consecuencia del consumo de este durante las horas del dfa). 
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Figura 23. Semivariograma experimental cruzado entre las variables clorofila "a" (54 datos) e IGC,(4). 
Informacion tomada en marzo de 1997 en la Cienaga Grande de Santa Marta. 

Una vez calculados los semivariogramas experimentales se ajusto el modelo lineal de 
corregionalizacion entre las variables CLORO 2 e IGC,(4) (tabla 7), el cual incluye efecto 
pepita puro y un modelo esferico. Con base en el modelo lineal de corregionalizacion se 
realizaron las predicciones, de la variable CLORO 2, a traves del metodo cokriging en 53 
sitios de muestreo (aquellos en los que fue eliminada inicialmente la informacion de clorofila 
"a") y se calculo la varianza de prediccion maxima, minima y promedio (tabla 8). Utilizando 
la informacion de los 54 datos de clorofila "a", se llevo a cabo prediccion en los 53 sitios 
restantes a traves del metodo kriging y se calculation nuevamente las varianzas de prediccion 
maxima, minima y promedio (tabla 8). Los resultados obtenidos por estos dos metodos 
fueron com parados con la varianza maxima de prediccion obtenida por Giraldo et al. (2001) 
para la valuable clorofila "a" utilizando la informacion completa (tabla 8). 


Tabla 7. Modelo de corregionalizacion ajustado a los semivariogramas experimentales (simples y cruzado) de 
las variables clorofila "a" (VI, 54 datos) e IGC(4) (V2, 114 datos). La informacion original fue medida en una 
muestreo realizado en marzo de 1997 en la C ienaga Grande de Santa Marta. 



Modelo Ajustado 

rM 

131.82 + 535.4 Esferico (8000) 


3.89 + 9.59 Esferico (8000) 

Yvlv2 (h) 

-1.18 - 18.70 Esferico (8000) 
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Tabla 8. Varianzas de prediction minima, maxima y promedio ((pg) 2 /l) para la variable clorofila "a", usando 
los metodos kriging y cokriging (con base en information auxiliar de la variable IGC(4)). Entre parentesis se 
encuentran la ganancia en precision respecto al metodo kriging con information incompleta. 


METODO 

VARIANZA DE 

PREDICCION 

MAXIMA 

VARIANZA DE 

PREDICCION 

MINIMA 

VARIANZA DE 

PREDICCION 

PROMEDIO 

Kriging con datos en 107 sitios de 
muestreo (information completa) 

379 (25%) 

0 


Kriging con datos en 54 sitios de 
muestreo (informacion reducida) 

506 (0%) 

0 

194.147(0%) 

Cokriging con datos en 54 sitios de 
muestreo para la variable clorofila 
"a" y 114 datos para la variable 
IGC/4) 

488 (4%) 

0 

190.06 (2.1%) 


De los resultados de la tabla anterior se concluye, como era de esperarse, que la 
reduccion en el numero de observaciones (reduccion del numero de sitios de muestreo) a 
cerca de la mitad, ocasiona un aumento en la varianza de prediccion de la variable 
considerada. No obstante dicho aumento es menor cuando se aplica el metodo cokriging 
utilizando como variable auxiliar - la variable IGC,(4). Al hacer la prediccion de la variable 
clorofila "a" con menos informacion empleando el metodo cokriging se gana un 4% en 
terminos de la varianza maxima y un 2.1 % , en terminos de la varianza promedio (promedio 
de 53 varianzas estimadas), de precision respecto al metodo kriging con datos incompletos. 

Teniendo en cuenta lo anterior y que la variable clorofila "a" presenta altos costos de 
muestreo (pasar de 107 sitios a 54 disminuye en mas del 200% (cerca de 2 millones de 
pesos) los costos (Giraldo et al., 2001) se recomienda en este caso la aplicacion del metodo 
de cokriging. Con la conclusion anterior no se indica que el numero optimo de puntos de 
muestreo para la variable clorofila "a" en el ecosistema de estudio debe ser 54, solo se 
muestra que en los casos en que se tiene informacion de variables auxiliaries como las aqur 
utilizadas es preferible el uso del predictor cokriging. Un estudio del numero optimo de 
sitios de muestreo requiere del diseno de una red de muestreo. Para lo anterior se puede 
consultar Giraldo et al (2001). 

5.5.2. Estudio Multivariado del Comportamiento Espacial de Variables Eisicoquunicas y 
Biologicas Medidas en el Estuario Cienaga Grande de Santa Marta. 

Se emplea la informacion de las variables salinidad, solidos en suspension (mg/1), 
nitritos (pmol/1), silicatos (pmol/lj y clorofila “a”(|ig/l), tomada de la base de datos descrita 
en las aplicaciones anteriores, para mostrar el uso de los componentes principales en la 
description de la distribution conjunta de estas. 

• Resultados y Discusion 

Con base en la matriz de correlacion de Pearson clasica (en distancia cero) y aplicando 
la metodologra descrita en la section 6.4.2. se generaron los ejes factoriales. Los dos 
primeros componentes explican el 71% de la variabilidad contenida en las 5 variables 
consideradas (tabla 9). Se consideran para el analisis solo estos dos componentes. La 
funcion de semivarianza cruzada entre los dos primeros componentes principales (Fig. 24 
revela que para distancias diferentes de cero hay asociacion entre los ejes principales y por 
consiguiente se puede concluir que no existe correlacion intrfnseca entre las variables 
originales. Lo anterior indica que un estudio exhaustivo de la correlacion entre las variables 
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requerirfa del calculo de matrices de corregionalizacion para diferentes distancias. En el 
documento se presenta solo el caso en que la distancia es cero, puesto que el proposito es 
identificar la distribucion espacial de las 5 variables de forma simultanea . 

Las variables de mayor importancia en la construccion del primer eje principal son salinidad 
(en sentido directo), silicatos y clorofila "a" (en sentido inverso). Analogamente las variables 
de mas peso en la construccion del componente dos son solidos en suspension y clorofila "a" 
en sentido directo y nitritos en sentido inverso (tabla 10). 

Tabla 9. Porcentajes de varianza explicados por los componentes principales generados con informacion de 
cinco variables fisicoqmmicas y biologicas medidas en el estuario Cienaga Grande de Santa Marta en marzo de 
1997. 


COMPONENTE 

VALOR PROPIO 

% DE VARIANZA 

% ACUMULADO 

1 

2.23 

44.610 

44.610 

2 

1.34 

26.962 

71.562 

3 

0.70 

14.124 

85.746 

4 

0.37 

7.426 

93.172 

5 

0.34 

6.828 

100.000 


Con la informacion georreferenciada de los dos componentes principales se realizaron 
analisis geoestadfsticos y se obtuvieron mapas de distribucion espacial (Figs. 25 y 26). El 
mapa de distribucion del componente uno (Fig. 25) indica que en gran parte del sistema el 
nivel de salinidad esta por debajo del promedio y que serfa esperable encontrar alii altas 
concentraciones de silicatos y clorofila "a". Hacia la zona norte del sistema se dan por el 
contrario magnitudes de salinidad por encima de su promedio y bajas concentraciones de 
silicatos y clorofila "a". En resumen podrfa pensarse, respecto a la informacion aportada por 
el primer eje principal, que la productividad biologica puede verse favorecida por altas 
concentraciones de nutrientes y bajos niveles de salinidad en sitios aledanos. 

La informacion aportada por el mapa de distribucion espacial del componente dos (Fig. 
26) confirma en gran medida lo descrito respecto a la distribucion espacial de la clorofila "a" 
dentro del sistema, es decir concentraciones altas en gran parte del cuerpo de agua y bajos 
niveles hacia la zona norte. En este caso se puede concluir que las magnitudes altas de esta 
variable estan asociadas a niveles altos solidos en suspension y , contrario a lo esperado, a 
bajos niveles en el ion nitrito. 

Los resultados descritos de forma conjunta a traves de la interpretation de las figuras 25 
y 26 son completamente coherentes con los reportados de forma univaiiada para las 
variables originales en la section 4.5.1. 

Tabla 10. Pesos de las variables en la construccion de los dos primeros componentes principales. Informacion 
original medida en marzo de 1997 en la Cienaga Grande de Santa Marta. 


VARIABLE 

COMPONENTE 1 

COMPONENTE 2 

Salinidad 

0.5060 

0.4150 

Solidos en Suspension 

-0.3468 

0.5668 

Nitritos 

0.2084 

-0.6219 

Silicatos 

-0.5049 

-0.3334 

Clorofila "a" 

-0.5703 

0.9140 
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Y (Ihl) 



Figura 24. Funcion de semivarianza cruzada entre los dos primeros componentes principales generados con 
information de algunas variables fisicoqufmicas y biologicas medidas en marzo de 1997 en el estuario Cienaga 
Grande de Santa Marta. 



Figura 25. Distribution espacial del primer componente principal generado con information de variables 
fisicoqufmicas y biologicas medidas en el estuario Cienaga Grande de Santa Marta en marzo de 1997. 
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Figura 26. Distribution espacial del segundo componente principal generado con information de variables 
fisicoqmmicas y biologicas medidas en el estuario Cienaga Grande de Santa Marta en marzo de 1997. 

5.5.3. Diseno de una Red de Muestreo para el Estuario Cienaga Grande de Santa Marta. 

En esta seccion se hace una aplicacion de la metodologfa propuesta por McBratney et al 
(1981) con el objetivo de disenar una red optima de muestreo para la CGSM (ampliamente 
descrita en secciones anteriores). En las dos ultimas decadas la CGSM ha venido dando 
muestras de deterioro y por ello se han implementado algunas obras civiles que buscan su 
recuperacion. Para el monitoreo de los cambios que se estan dando en el ecosistema se hace 
necesario contar con un conjunto de sitios de muestreo que haga posible lograr una vision 
integrada del comportamiento de las principales variables que rigen sus procesos de 
productividad. 

Se analizaron datos en la superficie de la columna de agua de las variables temperatura, 
salinidad, seston total, profundidad, silicatos, clorofila, oxfgeno disuelto, nitritos y clorofilas 
“a” y "c" tornados en los mismos puntos de muestreo descritos en la seccion 3.3. 

Se simularon redes de muestreo con cuadrfculas de 4, 9, 16, 25 y 36 km 2 , 
respectivamente (Fig. 27) y se estimaron las correspondientes varianzas de prediccion de 
cada variable en cada epoca tomando como base los modelos de correlacion espacial 
estimados en la seccion 3.3. 

La comparacion del error estandar de prediccion y de los costos asociados al muestreo 
de cada variable en cada red, posibilito el establecimiento de un conjunto de sitios de 
muestreo optimo bajo estos dos criterios. 
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Figura 27. Redes de muestreo bajo las cuales se hicieron las estimaciones de las varianzas de prediction de 
cada una de las variables consideradas, asumiendo los modelos de semivarianza estimados. Distancias entre 
puntos de muestreo: a) 2000 m; b) 3000 m; c) 4000 m; d) 5000 m y e) 6000 m. 

• Resultados y Discusion 

Si bien es posible que en la fecha del muestreo se estuviese dando un fenomeno de 
intervencion debido a los cambios climatologicos dados en el ano inmediatamente anterior a 
este, para los propositos del trabajo esto no resulta ser un inconveniente puesto que de hecho 
se asume que el establecimiento del conjunto optimo de puntos de muestreo no depende de 
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la magnitud de la variables sino de la estructura de dependencia espacial presente en la 
region de estudio para cada una de ellas. 

Como era de esperarse los errores estandar de prediccion estimados son crecientes en 
funcion de la distancia que existe entre los puntos de muestreo (tabla 11). La variable en la 
cual se consigue mayor ganancia relativa en precision al pasar de la red menos densa (Fig. 
27 e), a la mas densa (Fig. 27.a) es la salinidad. En dicha variable se consigue aumcntar en 
un 35 % la precision (tabla 12). Otras variables como temperatura, oxigeno disuelto, silicatos 
y clorofila “a” tienen aumentos en precision que oscilan entre 15.9 y 23.8 % (tabla 12). Por 
ultimo en las valuables profundidad, nitritos, seston total y clorofila "c", solo se consigue 
aumentar la precision en porcentajes que estan entre el 5,7 y el 10.1 % (tabla 12). 
Obviamente si se comparan las redes intermedias, redes con distancias entre puntos de 
muestreo de 3000, 4000 y 5000 m (Fig. 27.b, 27.c y 27.d), con la red de 6000 m (Fig. 27.e), 
resultan mucho menores los incrementos relativos en precision (tabla 12). 


Tabla 11. Error estandar (raiz cuadrada de la varianza) de prediccion maxima de cada variable considerada 

para redes de muestreo con cuadriculas de 4, 9, 16, 25 y 36 km 2 . 






TAMANO DE LA RED 


Variables 

(Distancia entre los puntos de muestreo en metros) 



2000 

3000 

4000 

5000 

6000 

Profundidad (m) 

0.2825 

0.2874 

0.2930 

0.3002 

0.307 

Temperatura (°C) 

0.6380 

0.6690 

0.7046 

0.7632 

0.8373 

Salinidad 

0.9096 

1.0511 

1.1676 

1.2965 

1.4075 

Oxigeno disuelto (mg/L) 

1.5145 

1.5917 

1.6752 

1.7977 

1.9431 

Seston total (mg/L) 

35.6363 

36.4021 

37.0459 

37.8076 

38.5197 

Nitritos (umol/L) 

0.2832 

0.2875 

0.2913 

0.2958 

0.3003 

Silicatos (umol/L) 

47.6524 

50.207 

52.3806 

54.6797 

56.6932 

Clorofila a (ug/L) 

19.4634 

21.2041 

22.5233 

23.5582 

24.2163 

Clorofila c (ug/L) 

6.1071 

6.2977 

6.4536 

6.6336 

6.7967 

Tabla 12. Ganancia en precision en porcentaje (cociente entre los respectivos 

errores estandar de prediccion) 

de cada una de las redes de muestreo (observada y simuladas) en 

cada variable respecto a la red de 6000 metros 

(la menos densa). 








TAMANO DE LA RED 


Variables 

(Distancia entre los puntos de muestreo en metros) 



2000 

3000 

4000 

5000 

6000 

Profundidad (m) 

8.0 

6.4 

4.6 

2.2 

0 

Temperatura (°C) 

23.8 

20.1 

15.8 

8.8 

0 

Salinidad 

35.4 

25.3 

17.0 

7.9 

0 

Oxigeno disuelto (mg/L) 

22.1 

18.1 

13.8 

7.5 

0 

Seston total (mg/L) 

7.5 

5.5 

3.8 

1.8 

0 

Nitritos (umol/L) 

5.7 

4.3 

3.0 

1.5 

0 

Silicatos (umol/L) 

15.9 

11.4 

7.6 

3.6 

0 

Clorofila a (ug/L) 

19.6 

12.4 

7.0 

2.7 

0 

Clorofila c (ug/L) 

10.1 

7.3 

5.0 

2.4 

0 
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Figura 28. Costos de muestreo de variables fisicoquimicas y biologicas en la Cienaga Grande de Santa Marta, 
segun diferentes espaciamientos entre sitios de muestreo (se asumen muestreos sistematicos de cuadrfculas). 

De otro lado si se estudian los costos de muestreo asociados a cada variable bajo cada 
una de las densidades muestrales (Fig. 28), se observa que existe considerable diferencia, 
con excepcion de las variables temperatura, profundidad y salinidad, entre la red de 2000 m 
y las restantes respecto a dichos costos. Para algunas de las variables (oxigeno disuelto, 
silicatos y clorofilas) pasar de la red de 3000 m a la 2000 m, implica incrementar el costo de 
muestreo de cada una de ellas en mas de $300000 

En conclusion para las variables temperatura y salinidad serfa mucho mas conveniente 
hacer un muestreo intensivo (red mas densa) dado que se consigue, comparando con la red 
menos densa, aumentar la eficiencia en porcentajes considerables (23 y 35%, 
respectivamente, tabla 12), con costos netos que se incrementan solo alrededor de $100000 
(Fig. 28). En la variable profundidad, si bien los costos de muestreo no se incrementan 
significativamente (Fig. 28), es mas aconsejable muestrear en la red menos densa dado que 
la eficiencia se incrementa en maximo un 7% al compararla con las restantes redes (tabla 
12). En las variables nitritos, seston total y clorofilas "a" y "c" hay poco aumento de la 
eficiencia al pasar de la red de 6000 m a otras con mayor numero de puntos (tabla 12) y por 
el contrario los costos, especialmente en la red de 2000 m, tienen aumentos considerables, 
lo que hace que se planteen las redes menos densas (5000 m y 6000 m entre puntos de 
muestreo) como las mas adecuadas para el seguimiento de estas variables. En las restantes 
variables (oxigeno disuelto, silicatos y clorofila “a”) es un poco mas compleja la decision 
dado que se obtienen aumentos considerables en los costos (Fig. 28), pero tambien 
incrementos de eficiencia (tabla 12). De todas formas es claro que se debe descartar en este 
caso la red con distancias entre puntos de muestreo de 2000 m dado que entre esta y la red de 
3000 m, la eficiencia relativa aumenta en un maximo del 8 % (tabla 12) con costos que se 
duplican o triplican para algunas variables (Fig. 28). 
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Apendice 


6.1. Indicador IGCi(P). 

Suponga que se tiene informacion dicotomica sobre P variables en n sitios de muestreo. 
La estructura de esta informacion se presenta a continuacion: 


Sitio 

X 

Y 

Vi 

v 2 ■ 


1 

x, 

Yi 

061 

061 ■ 

•• 061 

2 

x 2 

y 2 

061 

061 ■ 

•• 061 

3 

x 3 

y 3 

061 

061 ■ 

•• 061 

n 

X n 

y n 

061 

061 ■ 

•• 061 


donde X y Y representan las coiTespondientes coordenadas de ubicacion geografica (grados, 
planas o cartesianas), V/, . . V p , son P variables que indican la presencia (1) o ausencia (0) 
de la caracterrstica. Si las variables observadas son cuantitativas se pueden catcgorizar baja 
diversos criterios. Uno de los mas empleados es el de comparacion de cada valor con su 
correspondiente mediana muestral. 

La informacion referente a tales variables puede presentarse en una matriz de la siguiente 
estructura: 


\ 7,2 

Vn ^22 

\JJnl Vn2 


'Hip 

^2P 




Si la variable j en el sitio i cumple 

el criteriode dicotomizacion 
en otro caso 


Se define el indicador del numero de "exitos" por sitio como: 

P 

= , \/i,i = l,---,n 

j=i 

y el indicador del numero de "exitos" par - la variable j como: 

n 

Y J lhj=V. j ,VjJ = l,2,-,P 

i=i 
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Ademas sea: 


K , =Z^;,> donde S tj =| 


j Sin. =1 
%~P Sirijj =0 


El indicador IGC/(P) en el sitio /-esimo se calcula como: 


donde: 


IGq{p) = P{P -2) + //,. + K, - h{p, 


h{P,ri,) = \ 


( a 

v ' 


- 2P 


0 

m ~ b ,n P ) 


Si ?],.=() 

Sirj ,, =1 

Si m = 7j i% - 2 y =2,3,---,P 


a 0 =0 

con: / m = 1,2,-■■ 

a,„ = a,„ , + [nr + m ) 

m m—J V / 


y 


b 0 =-i 

b m = K-, + (m-l) ’ 


m — 1, 2,-■■ 


En la tabla 13 se presentan los valores los valores de a m , b m y h(p, ij t ) necesarios para el 
calculo de la variable IGQ(P) cuando el numero de presencias de especies esta entre 2 y 15. 

Tabla 13. Valores de a m y b m en la ecuacion 6 para rjj. = 2, 3, . . ., 15. P es el numero de variables de tipo 


presencia - ausencia consideradas. 


Pi. 

m 

d m 

K 

1ET 

3 

1 

<3- 

3 

2 

0 

0 

-1 

P 

3 

1 

2 

-1 

(2 + p) 

4 

2 

8 

0 

8 

5 

3 

20 

2 

(; 20-2 p ) 

6 

4 

40 

5 

0 40- 5 p ) 

7 

5 

70 

9 

(70 -9p) 

8 

6 

112 

14 

{ I12-14p ) 

9 

7 

168 

20 

(l68 — 20p) 

10 

8 

240 

27 

[240- 27p) 

11 

9 

330 

35 

[330-35p) 

12 

10 

440 

44 

[440- 44 p) 

13 

11 

572 

54 

[572- 54 p) 

14 

12 

728 

65 

[728-65p) 

15 

13 

910 

77 

[910-77 p) 
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Para ilustrar el calculo de la variable IGC(P), se presenta en la tabla 14 el caso particular de 
P = 4. Se utilizan todos los posibles arreglos de ceros y unos que se pueden formal - con 4 
datos. 


Tabla 14. Calculo de los valores de la variable IGQ(4) con todos los posibles arreglos formados con 4 
sfmbolos de dos tipos (ceros y unos). 


ARREGLO 

7,7 

7 a 

7 i3 

7 'a 

7 

S,i 

s a 

Si3 

s i4 

K, 

P(P -2) 

IGC(P), 

1 

0 

0 

0 

0 

0 

-4 

-4 

-4 

-4 

-16 

8 

0 

2 

1 

0 

0 

0 

1 

1 

-3 

-3 

-3 

-8 

8 

1 

3 

0 

1 

0 

0 

1 

-3 

2 

-3 

-3 

-7 

8 

2 

4 

0 

0 

1 

0 

1 

-3 

-3 

3 

-3 

-6 

8 

3 

5 

0 

0 

0 

1 

1 

-3 

-3 

-3 

4 

-5 

8 

4 

6 

1 

1 

0 

0 

2 

1 

2 

-2 

-2 

-1 

8 

5 

7 

1 

0 

1 

0 

2 

1 

-2 

3 

-2 

0 

8 

6 

8 

0 

1 

1 

0 

2 

-2 

2 

3 

-2 

1 

8 

7 

9 

1 

0 

0 

1 

2 

1 

-2 

-2 

4 

1 

8 

7 

10 

0 

1 

0 

1 

2 

-2 

2 

-2 

4 

2 

8 

8 

11 

0 

0 

1 

1 

2 

-2 

-2 

3 

4 

3 

8 

9 

12 

1 

1 

1 

0 

3 

1 

2 

3 

-1 

5 

8 

10 

13 

1 

1 

0 

1 

3 

1 

2 

-1 

4 

6 

8 

11 

14 

1 

0 

1 

1 

3 

1 

-1 

3 

4 

7 

8 

12 

15 

0 

1 

1 

1 

3 

-1 

2 

3 

4 

8 

8 

13 

16 

1 

1 

1 

1 

4 

1 

2 

3 

4 

10 

8 

14 


Se observa en la tabla anterior que el IGC,(4) es una variable discreta monotona 
creciente con valores entre 0 y 14. Para este caso los valores del indicador definen 
puntualmente, excepto cuando el IGC,(4) igual a 7, lo ocurrido respecto al numero de unos y 
a la posicion de los mismos dentro de las sucesiones. Valores del IGC,(4) entre 1 y 4 indican 
que hubo un solo uno en la sucesion y cada numero revela la posicion que este ocupa en la 
misma. Valores entre 5 y 9 corresponden a sucesiones en las que hubo dos unos, con IGC,(4) 
igual a 5 cuando los dos unos estan en las primeras dos posiciones de las sucesion ordenada 
y a 9 cuando estan en las dos ultimas. Los valores 6, 7 y 8 reflejan la transicion de las dos 
primeras a las dos ultimas posiciones. Valores entre 10 y 13 indican que hubo tres unos y 
cada uno de estos valores corresponde a una unica sucesion. El valor del IGC,(4) sera 10 
cuando los unos esten en las tres primeras posiciones e igual a 13 cuando esten en las tres 
ultimas. Los valores 0 y 14 se obtendran cuando en la sucesion no haya unos o todos los 
valores sean iguales a uno, respectivamente. 


6.2. Algebra de Matrices. 

La gran mayorfa de metodos estadfsticos, incluyendo la geoestadfstica, pueden ser 
tratados de forma mucho mas sencilla a traves del uso del algebra de matrices. Por esta razon 
es util, si no esencial, tener un cierto conocimiento mmimo de esta area de las matematicas. 
Lo anterior es cierto siempre y cuando el interes sea usar los metodos como una herramienta. 
La notacion del algebra matricial algunas veces puede resultar desanimante. Sin embargo, no 
es diffcil entender sus principios basicos. 
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6.2.1. Matriz 


Una matriz A de tamano (mxn) es un arreglo rectangular de m filas con n columnas. 

a ll a 12 a ln ^ 

a 21 a 22 a 2n 


V a ml a m2 • a mn' 


6.2.2. Suma y Producto de Matrices 

El procesos aritmetico de adicion, sustraccion, multiplicacion y division tiene sus contraparte 
con matrices. Si A y D son dos matrices de orden 3x2, entonces su suma se define como: 



r a ll 

a 12^ 


f d ll 

d 12^ 


( a ll + d ll 

a 12 + d 12^ 

A + D = 

a 21 

a 22 

+ 

d 21 

d 22 

= 

a 21 + d 21 

a 22 +d 22 


^ a 31 

a 32/ 


U32 

d 232 


^ a 31 + d 31 

a 32 + d 32/ 


En el caso de la multiplicacion se debe cumplir que el numero de columnas de la primera 
matriz sea igual ala numero de filas de la segunda. 

A Z a li b il Z a li b i2 Z a li b i3^ 



^ a ll 

a 12^ 

f 

A • B = 

a 21 

a 22 

• 

V 


^31 

a 32/ 



b ll 

b 21 


b 12 

b 22 


b 13^ 

b 232 


Z a 2i b il Z a 2i b i2 Z a 2i b i3 

Z a 3i b il Z a 3i b i2 Z a 3i b i32^ 


6.2.3. Inversa y Determinante de una Matriz. 

Si k es un numero, es cierto que k x k 1 = 1. De forma similar si A es una matriz 
cuadrada (numero de filas igual al numero de columnas) su inversa es A' 1 , donde AA" 1 = 
A*A = I, con I igual a la matriz identica (matriz de unos en la diagonal y cero por fuera de 
ella). Un ejemplo de matriz inversa es: 


"2 l) 

-1 

^ 2/3 -1/3' 

vl 22 


1-1/3 2 / 3 y 


Esto puede comprobarse observando que: 


"2 

vl 


2y 



f 2/3 

-1/3) 


fl 

0^ 

• | 







V-l/3 

2/3 j 


lo 

b 
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la inversa de una matriz 2x2, si existe, puede determinarse facilmente por medio del 
siguiente calculo: 


" a ll 

a 12^ 

-1 

f a 22/ 

/A 

—a 12/^1 
/A 

U 21 

a 22^ 


' a 2l/ 

a ll/ 


v /A /A J 


Donde A = ana 22 - ai 2 a 2 i. La cantidad A es llamada el determinants de la matriz. Claramente 
la inversa no esta definida si el determinante es igual a cero. Con matrices grandes el calculo 
de la inversa es tedioso y se debe usar un programa de computo para realizarlo. 

6.2.4. Valores y Vectores Propios. 

Dada una matriz A de orden (n x n), si existe un vector x (n x 1) y un numero X tal que 

Ax = Xx. 6 (A - ?d)x = 0 


donde I es la matriz identica de orden (n x n) y 0 es un vector (n x 1), entonces se llama a X 
y x, respectivamente, valor y vector propio de la matriz A. Pueden encontrarse hasta n 
valores propios y hay tantos vectores propios como valores propios se encuentren. Los 
valores de X deben satisfacer que el determinante de A - ?il = 0. Los vectores propios se 
calculan despues de reemplazar los valores propios encontrados en la expresion Ax = Xx. A1 
igual que con la inversa, para matrices grandes se debe emplear un software especializado 
para su obtencion. A continuacion, a manera de ilustracion, se realiza el calculo de los 
vectores y valores propios de una matriz de orden 2x2. 


Sea A = 


^6 

v3 


3^ 

4y 


, entonces 


|A - Xl\ = 0 : 



(6 33 | 

f 1 

03 

0^ 

-x\ 

v3 4 1 

lo 

lJ 


= 0 


% ^ 
v3 4y 


r X (C 

v0 Xy 


0 6-X ) 3 

3 (4 - X) 


0 


(6 - X)(4 - X) -9 = 0 
X 2 - 10X + 15 = 0 


x = 


-b ± 



- 4ac 


2a 
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2 


X = 8.1623, X = 1.8377 


Para cada valor propio existe un vector propio, el cual se obtiene reemplazando el valor 
propio correspondiente en la primera expresion de la pagina anterior y usando la condicion 
de que los respectivos vectores propios esten normalizados. 


Un vector x 


VX22 


se dice que esta normalizado si satisface que Jx j + x 2 = 1. 


Teniendo en cuenta lo anterior se calculan los vectores propios de la siguiente forma: 
(A - Xl)x = 0 

3 Yxi"! fO) 

" 1 Vx 2 2 


(4 -X)) 


loj 


(6-X)xi + 3x 2 = 0 
3x^ + (4 -A,)x 2 = 0 


Restando las dos ecuaciones anteriores y factorizando, obtenemos: 

x 1 (6-X-3) + x 2 (3-4 + X) = 0 
x 1 (3-X,) + x 2 (-1 + ^) = 0 



Entonces para X = 8.1623 y X = 1.8377 se tiene respectivamente: 

xi = 1.3847x2 y xi = -0.7207x2 . Ahora utilizando la restriccion de que los vectores 
esten normalizados se obtiene: 


xf =(l.3847) z l-xj 


X 2 + (1.3847) 2 x^ = (1.3847) 2 
Xj (l+ 1.3847 2 j = (1.3847) 2 


(1.3847) z 
' I o 

1 + 1.3 847 2 


X 1 


1.3847 


1 +1.3847" 


0.8107 


Reemplazando el valor de xi, obtenemos que x 2 


X 1 


0.8107 


1.3847 1.3847 


0.5855. 




Luego el vector propio asociado al valor propio X = 8.1623 es 


U22 


f 0.8107" 
0.5855y 


Efectuando un procedimiento similar se puede comprobar que el vector propio asociado al 


valor propio X = 1.8377 es 


'xj" 

U27 


-0.5847 
0.8113 ) 


En resumen dada la matriz del ejemplo entonces se puede comprobar que: 


'cn 

so 

1 _ 


( 8.1623 0 " 

YO.8107" 


^0" 

b 4j 


v 0 8.1623y 

vO.5855; 


la 


y, con el segundo valor y vector propio, que 


1 

o\ 


Y.8377 0 " 

Y-0.5847" 


(0) 

b \) 


v 0 1.8377y 

l 0.8113 y 


loj 


6.3. Conceptos de Probabilidad 

A continuacion se presenta una revision no exhaustiva y a manera introductoria de conceptos 
basicos de la teorfa de probabilidades. Un estudio profundo y formal de estos se puede hacer 
en Mood et al (1963) . 

6.3.1. Variable Aleatoria 

Si X es una funcion que le asigna a cada uno de los resultados de un experimento aleatorio 
(aquel cuya respuesta no puede ser establecida de antemano) un numero real, entonces X se 
llama una Variable Aleatoria. Estas pueden ser discretas o continuas. 

6.3.2. Funcion de Probabilidad 

Si X es una variable aleatoria discreta. Se llamara a f(x) = P (X = x) funcion de probabilidad 
de la variable aleatoria X, si satisface las siguientes propiedades: 

i. f(x)>0 \/xeR x 

11. £/(*)=;. 


Si existe una funcion f(x) tal que: 


i. f{x) >0, — °° < x <°° 

ii. | f(x)dx = l 

iii. p{a < X b)= f b f(x)dx para cualquier a y b, entonces f(x) es la funcion de densidad de 

J Cl 

probabilidad de la variable aleatoria continua X. 
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La funcion de probabilidad acumulada, notada como F(x), es igual a p(x < x) y se evalua a 
traves de una sumatoria o de una integral dependiendo de si X es discreta o continua. 


6.3.2.I. Valor Esperado y Varianza 

Si X es una variable aleatoria, el valor esperado de una funcion de la variable aleatoria X, 
g(x)esta dado por: 


E{g{x)) 


XsM/M 

f g{x)f{x)dx 

J —oo 


X discreta 
X continua 


como caso particular, 


E{x) = fi 


X 

[ xf(x)dx 

J —oo 


X discreta 
X continua 


La varianza de la variable aleatoria X esta definida como: 


V(X ) = a 2 


E{X-juf 


Yj( x ~v ) 2 fi x ) 

j {x-juf f{x)dx 


X discreta 
X continua 


La rarz cuadrada de la varianza se denomina desviacion estandar y se denota por <J . 

Se cumple que: 

1 . E(aX ) = ciE(x ) , con a constante 

2. E(aX +b) = aE(x ) + b, con ay b constantes 

3. v{aX) = a 2 v{x ) y a constante 

4. y(x) = £(x 2 )-[£(x)] 2 

6.3.2.2. Funcion de Probabilidad Binomial y Normal. 

Modelo Binomial 

Suponga que hay un experimento que consiste en examinar n individuos y evaluar o medir 
en cada uno de ellos si tienen o no una caracterfstica dada (solo hay dos posibles 
resultados).Sea p la probabilidad de "exito" y q = 1-p la de "ffacaso" en cada uno de los n 
ensayos. Se asume que esta probabilidad es constante en cada uno de ellos. 
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Sea X = Numero de exitos en los n ensayos, entonces asumiendo conocido p entonces es 
posible establecer las probabilidades de ocurrencia de cada evento mediante la siguiente 
ecuacion, denominada modelo de probabilidad binomial. 


P(X=x) = 


P x 0 - P) 


x = 0,l,2,...,n 




En este modelo: 


p = E(X) = np 

a 1 = V(X) = np{l-p) 


Modelo Normal 


El modelo de probabilidad normal (Gaussiano) es util para encontrar las probabilidades 
asociadas a eventos de variables aleatorias cuyas distribuciones de frecuencias son simetricas 
alrededor del valor promedio. Algunos ejemplos de este tipo de variables aleatorias son los 
siguientes: 

Sea p el valor promedio de la variable (E(X)) y o 2 su correspondiente varianza (V(X)), entonces 
las probabilidades de ocurrencia de eventos asociados a los posibles resultados de la variable 
estudiada pueden ser encontrados usando la siguiente expresion, llamada modelo de 
probabilidad normal 

b , 1 2 

P(a < X < b) = f ._ e 1 ’ dx- 

{ 42na 


Obviamente resultana muy dispendioso tener que calcular estas integrales para cada valor de 
a, b, p y o . Por esta razon se acude a un procedimiento llamado estandarizacion, el cual 
consiste en hacer la transformacion z _ x ~JL ■ La variable anterior tendra (si la distribucion 

O 


de frecuencias de X se ajusta a un modelo de probabilidad normal con media p y varianza 
O') una distribucion de frecuencias que se ajusta a un modelo de probabilidad normal con 
media cero y varianza uno, es decir que: 


P( a< X <b) = 


a- fi 


<Z< 


b - p 


^■1 

= (4; <Z<Z 2 )= J 


1 “V Z 


V2T 


La ecuacion anterior tambien puede resultar dificil de evaluar, sin embargo para cualquier 
valor de a, b, p y o las correspondientes probabilidades pueden hallarse, sin necesidad de 
resolver la integral, empleando la tabla de distribucion acumulada normal estdndar que 
aparece en los textos de estadfstica. 


6.3.3. Funcion de Probabilidad Bivariada. 


Si X y Y son dos variables aleatorias discretas. La probabilidad deX = xyY = y esta 
determinada por la funcion de probabilidad bivariada f(x,y) = p[x = x, Y = y] donde : 


i. f(x,y)>0,\/x,yeR x ,R Y 
ii- ^Yjf{x,y) = l 

x y 
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Si existe una funcion f(x,y) tal que la probabilidad conjunta: 

p\a<X <b,c<Y <d]=\ f f(x, y)dydx 

Ja Jc 

para cualquier valor de a, b, c y d en donde f(x,y)>0, -oa <x ,y<°° y J°° J°° f(x,y)dydx = l , 
entonces f{x,y) es la funcion de probabilidad bivariada de X y Y. 

La funcion de probabilidad acumulada f(x, y) es igual a p[x < x,Y < v] y se evalua a traves de 
una doble sumatoria o de una doble integral dependiendo de si las variables aleatorias son 
discretas o continuas, respectivamente. 

6.3.3.I. Funcion de Probabilidad Marginal 

Si X y Y son dos variables aleatorias con funcion de probabilidad conjunta f{x,y). Las 
funciones de probabilidad marginales de Y y Y estan dadas por 


/(*) = X f( x ’ y ) 

y 

f(y)=Yjfi x ’y) 


o por 


fi x )=\_ aa f{ x >yMy 

f(y)=ljUy dx ) 


si X y Y son variables aleatorias discretas 


si X y Y son variables aleatorias continuas 


6.3.3.2. Funcion de Probabilidad Condicional 


Sean X y Y dos variables aleatorias con funcion de densidad conjunta f{x,y). La funcion de 
probabilidad condicional de la variable aleatoria X, denotada por f(x/ y), para un valor fijo y 
de Y, esta definida por: 


fix/y) 


fjx.y) 

fiy) 


, donde fiy) es la funcion de probabilidad marginal de Y de manera tal que 


fiy)>o . 


De manera analoga, la funcion de probabilidad condicional de Y para un valor fijo x de X se 
define como: 
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f{y/x) = 
f(x)>0 . 


/U y) 

/(■*) 


, donde f(x) es la funcion de probabilidad marginal de X de manera tal que 


6.3.3.3. Independencia Estadistica. 

Sean X y Y dos variables aleatorias con funcion de densidad conjunta f{x,y). X y Y son 
independientes si y solo si: 

f(x,y)=f(x)f(y) 

donde f(x) y f(y) son las funciones de probabilidad marginales. 


6.3.3.4. Valor Esperado, Varianza y Covarianza 

Sean X y Y dos variables aleatorias que se distribuyen conjuntamente. El valor esperado de 
una funcion de X y Y, g(x, y), se define como: 


E{g{X,Y)) 


Y.Hs{x,y)f{x,y) 

Loo Loo g( x > y )f( x ’ y ) dydx 


si X yY son discretas 
si X yY son continuas 


La covarianza entre X y Y, denotada por Cov (X, Y), se define como: 

E[{x - fix ){Y fly )] = E{XY - Xfly - Yfl x + fl x JU Y ) = e{xy) - E{x) e{y) 

donde Mx y My representan los valores esperados de X y Y respectivamente. 

Si la covarianza de X y Y se divide por el producto de las desviaciones estandar de X y Y, 
el resultado es una cantidad sin dimensiones que recibe el nombre de coeficiente de 
correlacion y se denota por p{x,Y ). 

„(x.r)= c °^- r ) 

a X G Y 


6.3.3.4.I. Propiedades del Valor Esperado y la Varianza. 


Si X y Y son dos variables aleatorias con densidad conjunta, entonces se cumple que: 


1 . e(x +y) = e(x) + e(y) 

2. v{x ±y) = v(x) + v(y)±2Cov{x,y) 
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3. V Y j a l X l = '£Y i a i a j Cov(x i ,Xj). 

\i=l y i=l j=l 

Observacion: Cov(x i ,X j )=Cov(x j ,X i ) y Cov(x i ,x i ) = v{x j ) 

Como caso particular: 

V{a,X, ± a 2 X 2 )= aMx,) + «fv(X 2 )± 2Cov(X l ,X 2 ) 

3. Si £(x) = £(y), entonces ^-E^x- y) 2 ]=N{x)+^v(y)-CoyIx,y). 

6.4. Algunos Metodos Estadisticos. 

6.4.1. Regresion Simple 

En el modelo de regresion simple se establece una relacion lineal entre la esperanza 
condicional de una variable aleatoria Y dados unos valores fijos de una variable X. 

Modelo Poblacional 

Yj = /3 0 + fit Xj + £j 

E{Y / X i ) = Y i = P 0 + fijXf 

Y;: i-esimo valor de la variable respuesta o dependiente en la poblacion 

V, :i-esimo valor de la variable predictora o independiente en la poblacion 

/3 0 y p, son parametros poblacionales que representan el intercepto y la pendiente, 

respectivamente 

£i : i-esimo error aleatorio en la poblacion. 

Supuestos del Modelo. 

1. E( £i ) = 0 

2. V( £j )=C7 2 

3. Cov(£;,£j) = 0 

4. e , ~ Nip, a 2 ) 

Modelo Muestral 

y t =Po + Pi x i + c 

yt = y, + c- 

y;: i-esimo valor de la variable respuesta en la muestra 
V, : i-esimo valor de la variable predictora . 
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Po y Pi son las estimaciones de los parametros con base en la informacion muestral. 
e x : i-esimo erro muestral. 

Estimacion de p 0 y (3, 

Uno de los metodos de estimacion de los parametros es el de mmimos cuadrados, que 
consiste en encontrar los estimadores que hacen minima la suma de cuadrados de los errores, 

es decir aquellos valores que hacen mas pequena ^ e? = X(ft - ft) . 

i=l i-1 


n n i \2 n 

Yj £i =X (ft - ft) = X (ft - Po - Pi x i Y • Derivando e igualando a cero se obtiene: 


i=l i-1 

n 

i=l 


i-1 




M =-2Y J {Y i -/3 0 -p I x i )=0 y -&— = -2£x i {Y i -p 0 -p 1 x i )= 0. 

°Po i=l °Pl i=l 


A1 simplificar las dos ecuaciones anteriores y distribuir las sumas se tiene: 


Z Y i = n fio + Pi Z 


i=l 


i=l 


Z x i Y i =aZa +AZ- 


i=7 


i=l 


i=l 


Las dos ecuaciones anteriores se conocen como ecuaciones normales. Dadas las 
realizaciones y 1 ,y 2j ■■■,y„ las ecuaciones pueden resolverse para encontrar los estimados de 
los parametros: 

Z -b = nfa + AZ*/ 

i=l i=l 

y = P 0 + Pix 
Po = y-PP 


Z x i =(y-Pi x \ Z x ‘ + A Za 2 


i=l 


i=l i=l 


Zal = 

(=7 


z* 

1=1 _ p I 1=1 


f n \\ 

Z-b 


n 

v y 


Za + AZ-b 2 

i=7 7=7 
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n n 


n 2> *X*«- 

T^i =-—-— fa 


/ \ 2 
n 


i=l 


i=l 


n n 


*Z*< 

Z i=7 i=l 
x iyi - 

A=—— 


Zta ~~ - v X- y) 


i=l 


5>f- 

i-1 


/ \2 
n 

5> 

V /=•/ / 

n 


' Eh -*) 1 


i=l 


Se puede demostrar que los errores estandar estimados de los estimadores de los parametros 
corresponden a: 




6.4.2. Analisis de Componentes Principales. 

El analisis de componentes principales es uno de los metodos multivariados mas 
difundidos, que permite la estructuracion de un conjunto de datos de multiples variables de 
una poblacion, cuya distribucion de probabilidades no necesita ser conocida (Lebart et al., 
1995). 

Se trata de una tecnica matematica que no requiere un modelo estadfstico para explicar 
la estructura probabilfstica de los errores. Sin embargo, si es posible suponer que la 
poblacion muestreada tiene una distribucion conjunta normal multivariada, podra estudiarse 
la significacion estadfstica de los componentes y sera posible utilizar la muestra 
efectivamente observada para efectuar pruebas de hipotesis, que contribuyan a conocer la 
estructura de la poblacion original, con un cierto grado de confiabilidad, fijado a priori o a 
posteriori (Pla, 1986). 

Los objetivos mas importantes del analisis de componentes principales son: 

i. Generar nuevas variables que puedan expresar la informacion contenida en el grupo 
original de datos. 

ii. Reducir la dimensionalidad del problema que se esta estudiando, como paso previo para 
futuros analisis. 

iii. Eliminar, cuando sea posible, algunas variables originales, en el caso de que aporten poca 
informacion. 

Este analisis se basa en una transformacion lineal de las observaciones originales. Esta 
transformacion es conocida en el campo del algebra vectorial como generacion de vectores y 
valores propios. Las nuevas variables generadas se llaman componentes principales y poseen 
algunas caracterfsticas estadfsticas deseables, tales como la independencia (cuando se asume 
la multinormalidad) y en todos los casos la no correlacion. Esto significa que si las variables 
originales no estan correlacionadas, los componentes principales no ofrecen ventaja alguna. 
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Generacion de los Componentes Principales 

Se ha dicho que los componentes principales tienen ciertas caracterfsticas que son 
deseables: 

a) Los componentes principales no estdn correlacionados. 

b) Cada componente principal sintetiza la maxima variabilidad residual contenida en los 
datos. Es decir, el primer componente sintetiza la maxima variabilidad posible en el 
conjunto de datos originales; el segundo componente sintetiza la maxima variabilidad 
restante, sujeta a la condicion de no correlacion con el primer componente, y asf hasta el 
p-esimo componentes. 

c) Cada componente contiene informacion de todas las variables pero en diferentes 
proporciones. 

Matricialmente se expresa la generacion de los componentes a traves de: 

y _ y t r\—1 

(nxp) ~ ^ (nxp)*~ / {pxp)^'(pxp) 

donde : 

Y: Matriz cuyas columnas representan las nuevas variables (componentes principales). 
Estas tienen la propiedad de ser no correlacionadas. 

X: Matriz de datos originales 

L: Matriz de vectores propios de: a) X T X, si X es la matriz de datos originales; b) S (matriz 
de varianzas y covarianzas) si X es centrada; c) R (matriz de correlacion) si X esta 
estandarizada. 

D: Matriz diagonal con valores en la diagonal iguales a la raiz cuadrada de los valores 
propios de X T X , S o R. 

La transformacion lineal para generar los componentes principales (matriz Y) se fundamenta 
en el proceso de diagonalizacion de una matriz, X T X , S o R., segun el caso, a traves del 
teorema de descomposicion del valor singular 
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